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Chapter 1 


Introduction. 


In the introductory section the author would like to convince the Reader 
that nonlinear chemical systems are interesting and worth studying, because 
they play a very important role in numerous phenomena occuring in living 
and nonliving nature. Moreover, this chapter contains examples of appli- 
cations of nonlinear chemical systems to information processing of various 
types. Further, the basic facts about chemical excitable systems are given 
and possible ways of experimental realizations of such systems are described. 
This information enables to specify the aims of this work and determines the 


background and motivation for the research described in Chapter 2. 


1.1 Why are nonlinear chemical systems im- 
portant? 


Nonlinear chemical reactions are of vital importance for all living creatures as 
well as for numerous events occurring in the nonliving world. If the reactions 
proceed in far from equilibrium conditions, they may lead to a variety of 


phenomena (for example oscillations, trigger waves, excitable pulses etc.), 
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which contribute to and reveal in processes like morphogenesis ([1], p. 593), 
formation of static and dynamical spatial patterns ([1], p. 327-590) and 
the cycles, observed on all levels of organization of the living matter, from 
single cells to whole organisms ([1|, p. 140). In particular, the existence of 
cycles can be explained by the oscillatory character of the underlying chemical 
reactions [1]-[4]. Together with the cyclic external perturbations (for example 
by the day/night cycle) they give grounds to so called ” biological clocks” [1], 
ie. the mechanisms which control rhythms of different activities of the living 
creatures. Thus, the biological and biochemical systems provide a ” natural” 
environment for nonlinear chemical reactions. Lots of examples can be found 
in any textbook on biochemistry (eg. [5]) or review reports on oscillations 
in chemical reactions (eg. [6, 7]). The majority of the biochemical cycles are 
very complicated, as there are thousands of complex reactions proceeding 
inside the living cells, most of which are catalyzed by enzymes and contain 
autocatalytic/nonlinear steps. Still, they provide a very attractive subject of 
research, because of their fundamental importance. 

Oscillations have also been discovered in chemical reactions proceeding 
outside living organisms. Historically first well described example is known 
as the Bray reaction [8] (catalytic decomposition of HO», during which the 
rate of Og production oscillates). In early 1950s periodical changes of con- 
centrations of reagents were observed by B. P. Belousov [9] in the reaction 
of oxidation of citric acid with bromate, catalyzed with cerium. This ob- 
servation surprised the community of ” orthodox” chemists. Further studies 
confirmed that the oscillations were really the result of the chemical reaction. 
Moreover, other examples of catalytic oxidation of organic substances with 
bromate were found [10, 11]. This group of reactions, called the Belousov 


- Zhabotinsky (BZ) reactions, became probably the most famous and most 
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intensively studied type of oscillatory chemical reactions in the last 30 years. 
In the meantime, oscillations of concentrations of reagents were found in 
numerous catalyzed or non-catalyzed reactions (the references [2, 3, 4] give 
many examples of such systems). 

All these reactions have attracted a lot of scientific interest and there are 
at least two reasons for that. First, these reactions are usually far simpler 
than enzymatic biochemical processes and yet qualitatively similar to them, 
as they reveal oscillations, excitability etc. Consequently, their cognition and 
description is much easier, but it gives grounds for better understanding of 
complex mechanisms of biochemical cycles. Second, in systems where chem- 
ical reactions appear together with the process of diffusion, a diversity of 
phenomena have been observed, which are interesting themselves, as a sep- 
arate subject of studies. They include: chemical waves [10, 11], travelling 
fronts and excitable pulses (examples further in this work), Turing structures 
[12, 13, 14], pattern formation [1], [15]-[31], synchronization [32, 33, 34] and 
resonances [35, 36, 37], to name only a few. The patterns and structures 
appearing in such systems are not only beautiful [15], but they also influ- 
ence processes of practical use, for example chemical reactions catalyzed by 
surfaces [38]-[41]. 

Experimental investigations of these phenomena are (in most cases) very 
complicated. To make them easier and - first of all - to understand them 
better, their proper theoretical (and mathematical) description is necessary. 
Thus, much effort was devoted to finding the models of the known nonlinear 
reactions and processes. As a remarkable part of my thesis concerns the phe- 
nomena occuring in systems based on the Belousov - Zhabotinsky reaction, it 
is worth mentioning that a commonly accepted mechanism of this particular 


type of reaction was proposed in 1972 by R. J. Field, E. Kórós and R. M. 
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Noyes [42]. It is known as the FKN mechanism of the BZ reaction. The FKN 
mechanism gave grounds for mathematical description of the system. Two 
simplified mathematical models, which are derived on the basis of the FKN 
mechanism, are the Oregonator model [43] and the Rovinsky - Zhabotinsky 
model [44, 45, 46]. They are very popular and frequently used in numerical 
studies of the BZ reaction. 

The results obtained within mathematical models are characterized by a 
high level of generality. They can be applied to any systems, to which a given 
model is applicable. In particular, the oscillatory or excitable systems appear 
in a variety of fields apart from chemistry. Investigations of thermochemical, 
electrochemical or biochemical systems [47], pulses in heart tissues [48]-[52], 
optical systems [53], or even ecosystems (population growth, [1], p. 1-107, 
disease spread |1|, p. 610) or economics are only a few examples. 

On the other hand, models invented in other fields were adapted by 
chemists. A good example is the well known FitzHugh - Nagumo model 
(FH-N), originally developed to describe the propagation of pulses of exci- 
tation in nerve tissues [54, 55], as a simplification of the Hodgkin - Huxley 
model [56]. The FH-N model has been used for investigations of chemical 
systems for many years [57, 58]. 

Interesting applications of oscillatory or excitable chemical reactions in- 
clude information processing (this area of research is briefly introduced in 
Section 1.3) or self-oscillating polymer gels [59]-|63], which may be used to 
construct artificial muscles. 

Undoubtedly, the investigations of nonlinear chemical systems have be- 
come a fast developing field of science, with interdisciplinary character. In 
this field theoretical studies and computer simulations are equally important 


as experimental works and both approaches are complementary. 
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1.2 Excitable systems. 


This section contains basic information about chemical excitable systems. It 
includes the definitions of excitability, refractory time and explains how trav- 
elling excitable pulses may appear in spatially distributed systems. Further 
the response of an excitable system to periodic perturbations is discussed, 
using a very general and simple model. It is shown that the dependency of 
the frequency of excitations of the system on the parameters of perturbations 
may have a devil’s-staircase-like form, due to non-zero refractory time of the 


excitable medium. 


1.2.1 What is an excitable system? 


A chemical system is called excitable if it has a stable stationary state of 
a special kind. The amplitude of a small perturbation of this state decays 
exponentially in time. However, if the amplitude of a perturbation exceeds 
a certain critical value, then during the evolution towards the stable state, 
the distance on the phase space between the time dependent state of the 
system and the stable state may become a few orders of magnitude larger 
than the distance between the stable state and the initial perturbation. Thus, 
perturbations with sufficiently large amplitude are amplified and the system 
returns to its steady state after much longer time than in the case of a 
small perturbation. The evolution after such perturbation forms a closed 
trajectory (so called excitation cycle) in the phase space. For a given type of 
perturbation one can introduce Smin as the minimum strength of perturbation 
which leads to systems excitation. 

The excitable systems have another important and interesting feature: 


just after an excitation they become refractory with respect to consecutive 
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perturbations. The minimal time after which the system can be re-excited 
is called refractory time (t,). Of course, t depends on the strength of the 
system’s perturbation s. 

As far as the mathematical description of such systems is concerned, 
excitability means that the nullclines of equations describing the system, as 
well as the vector fields driving the variables, are qualitatively similar to those 
presented in Figures 1 and 2 (Section 2.1; the figures contain also examples of 
excitation cycles). Of course, any translations, rotations and mirror images 
of the presented nullclines together with the vector field are allowed. What 
is important is that the nullclines have only one common point, located near 
the extreme of the S-shaped curve, as well as the orientation of the vector 
field. If the nullclines intersect in three or more points (a system has two or 
more stable stationary states), then the excitable (non-monotonic) switching 
between the stable states may be observed ({1], p. 335). To describe an 
excitable system a set of minimum two autonomous differential equations 
(two-variable system) is necessary. 

If the system is distributed in space and locally perturbed, then a local 
excitation may propagate (by diffusion) to the neighboring areas, while the 
points initially excited relax to the stationary state. Thus, for the proper 
choice of parameters travelling pulses may appear in the system (one pulse 
per one supercritical perturbation). The pulses are stable. After initiation 
they develop into stationary forms, which do not depend on the method of 
initiation. 

The refractory time of the medium limits from above the frequency at 
which excitable pulses may appear at a certain point of a spatially distributed 
system. Once a pulse has passed through a part of the system, then at least 


the time t is necessary before another pulse would be able to propagate 
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in this region. If the next pulse arrives earlier it cannot propagate and 
disappears. The situation here is similar to burning the grass: once the grass 


has been burnt it needs some time to grow again before it can be burnt. 


1.2.2 Devil’s-staircase-like response in a” naive” model 


of a periodically perturbed excitable system. 


Let us consider an excitable system which is periodically perturbed. Observ- 
ing consecutive excitations of the system and comparing their frequency with 
the frequency of perturbations one may see that the frequency of excitations 
is equal to the frequency of perturbations, or lower. The difference between 
the frequencies of perturbations and of successful excitations can emerge in 
every excitable system, just because of the non-zero refractory time. If the 
time shift between consecutive perturbations is shorter than the refractory 
time of the perturbed medium (for a given strength of perturbations), then 
some of the perturbations have to occur at moments, when the system is in 
the refractory state and cannot be excited. Thus, for a high frequency of 
perturbations, not all of them lead to excitation of the system. 

Let us consider an excitable system characterized by the refractory time 
t,, schematically shown in Figure 0A. The horizontal axis of Fig. 0A presents 
the strength of the perturbation, given in Spin units (Smin is the minimum 
strength of perturbation which leads to system’s excitation). The vertical 
axis gives the refractory time (in arbitrary units). Let us take into account 
two successive perturbations, occurring at times 0 and t. Let us fix the 
strength s of these perturbations. At time 0 every perturbation charac- 
terized by s (s/Smin>1) excites the system. A successful excitation by the 
second perturbation depends on the values of t and s. The line t(s) in Fig. 


OA corresponds to the minimal time shift after which the second perturba- 
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Time between perturbations t (arb. units) 


Strength of perturbation s (in Smin units) 


Figure 0A. The ”naive” model of a perturbed excitable system. The proba- 
bility that the second perturbation excites the system is shown as a function 
of the strength of perturbation s in Spin units (x - axis) and the time sepa- 
rating perturbations t (arbitrary units, y - axis). The set of parameters for 
which the second perturbation excites the system is located above the line. 
For a given s the line separating excited by the second perturbation from 
unexcited ones shows t,(s). Note that for times t<t, and for s/Smin<1 the 


system cannot be excited. 
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tion of the strength s is able to excite the system. Thus, if the perturbation 
is strong enough, then re-excitation occurs when t>t,. The shape of the 
t,(s) dependency shown in Fig. OA seems quite natural, for the following 
reasons: First, if s is too small (s/smin<1), the system is not excited by the 
second perturbation, nor by the first one, no matter how large is the time 
shift between them. Second, after a successful excitation the system relaxes 
towards the steady state. It can be excited again if it has reached the ap- 
propriate neighborhood of the steady state before the second perturbation 
occurs. This neighborhood is larger (and reached faster) for stronger pertur- 
bations. Consequently, for s/smin>1 one should expect that t, gets smaller 
when perturbation’s strength increases. Third, even if the perturbations are 
very strong, the system cannot be excited if it is too far from the steady state. 
As a consequence, the curve t,(s) becomes less and less steep for larger s. 

For a simplified analysis let us assume that a state of the system is not sig- 
nificantly changed after an unsuccessful perturbation and the probability of 
excitation shown in Fig. 0A remains valid if a periodic sequence of perturba- 
tions is considered. Now ¢ denotes the time from the last perturbation which 
excited the system. Let us consider a periodic sequence of perturbations at 
times k-t,, kEN, characterized by a fixed strength s. It is easy to see that if 
tp>t,(s), then each perturbation excites the system. If t,(s)/2<t,<t, (s) then 
every second of them excites the system, if t.(s)/3<t,<t,(s)/2, then every 
third perturbation does it and so on. Plotting the firing number f (ie. the 
ratio between the frequencies of excitations and perturbations) as a function 
of t (Figure 0B), one obtains a plot similar to devil's staircase [64]. Another 
staircase-like plot can be obtained by representing the firing number as a 
function of s for a fixed t, (cf. Figure 0C). 


Of course, the ”naive” model that has just been discussed is extremely 
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Figure OB. The firing number as the function of the time interval be- 
tween consecutive perturbations of the system (in t,/t,(s) units), for a fixed 
strength of perturbation s. This result has been obtained for the ” naive” 


model of a perturbed excitable system. 
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Figure OC. The firing number as the function of the strength of perturbation 
s, for a fixed tp. This result has been obtained for the ”naive” model of a 


perturbed excitable system. 
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crude, but it predicts that devil’s-staircase-like dependence of firing num- 
ber on the period of perturbations should be widely observed in excitable 
systems. It is true indeed. Such dependence can be found in a number of ex- 
periments concerned with periodic perturbations of a homogeneous excitable 
system, the results of which were published by M. Marek and collabora- 
tors [65|-[68]. Another staircase-like dependence of the firing number was 
observed by A. Toth, V. Gaspar and K. Showalter in spatially distributed 
excitable systems linked with a capillary tube [69]. The excitable waves were 
generated in one system and they might excite the other one via reagents’ 
transport in the capillary. In this experiment the diameter of the capillary 
controlled the strength of perturbation. If it was small, the refractory time of 
the perturbed system was much longer than the refractory time in the system 
where waves were originally generated. As the result only a fraction of gen- 
erated waves excited the system after passing through the capillary and the 
firing number as the function of entering wave period had a devil’s-staircase- 
like dependence (cf. Fig. 6 in [69]). Similar phenomenon was observed by 
K. Suzuki, T. Yoshinobu and H. Iwasaki [70] in the ferroin catalyzed BZ 
system. In their experiment a train of excitable pulses, travelling on a fer- 
roin covered membrane, arrived at a gap (a narrow stripe of area without 
ferroin). The frequency of pulses on the other side of the gap depended on 
the gap’s width and was a certain fraction of the frequency of the arriving 
train of pulses. Moreover, recent studies indicate that it is possible to excite 
a spatially distributed excitable system by a local periodic modulation of one 
of its parameters [71]. Although all the perturbations are subcritical (ie. a 
single perturbation of that kind never leads to an excitation cycle), an ap- 
propriately tuned frequency of periodic perturbations goes in resonance with 


the damped oscillations around the stable fixed point and leads to regular 
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excitations of the system [71]. Also in this case the firing number (the ratio: 
number of excitations to number of perturbations in the same time inter- 
val) as a function of the modulation frequency reveals a devil’s-staircase-like 


structure. 
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1.3 Information processing by chemical sys- 
tems. 


The idea of information processing with purely chemical devices is attrac- 
tive for both scientific community and leading hi-tech companies. Two most 
promising ways of constructing such devices come from the molecular chem- 
istry approach (the information is processed by individual molecules) [72] 
and from the nonlinear chemical reactions. 

Logical functions and other kinds of information processing tasks were 
implemented in homogeneous (stirred) reactors, in which a reaction possess- 
ing two or more stable stationary states proceeds. Once a logical value is 
assigned to all of these states, the reactors can be applied to processing in- 
formation [73]. A good example is the idea of Okamoto and collaborators 
[74], who proposed to use a cyclic enzyme system, in which two reactions are 
coupled by a shared reagent. Such system can be controlled (or stimulated) 
by the supply of the substrate precursors of the two main reactions (consid- 
ered as the two *” input signals”). The level of the shared reagent, which has 
two stable stationary states, corresponds to the ”output signal”. This sys- 
tem behaves like a neural cell and because of that it is called a ” biochemical 
neuron”. The homogeneous reactors coupled by diffusion can be arranged 
in networks, and used eg. for pattern recognition [75] or building multilayer 
biochemical neuron networks, which may perform tasks of a high complexity 
(74, 76]. 

Two dimensional photosensitive BZ systems are also suitable for direct 
image processing. The BZ system should be prepared in the oscillatory 
regime. In this regime the color of the BZ solution depends on the point of the 


limit cycle and changes periodically in time. The period of oscillations can 
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be influenced by external illumination. In particular, if an image is projected 
onto such system, then a chemical ” representation” of the image is produced 
in the system, as the period changes from point to point, according to the 
light intensity. Once the light is switched off, different areas of the system 
start to evolve towards the new limit cycle with different initial conditions. 
Because of that the initial variances in illumination of neighboring areas 
may be increased (contrast enhancement) or decreased (contrast reduction) 
during the evolution. When the phase of the limit cycle changes by m, then 
the negative image is obtained. Alternatively, if a high intensity of light is 
used, then in some regions (the brightest parts of the image) the system 
may be switched from the oscillatory to excitable regime. This regime may 
be applied to determining or smoothing contours of the originally projected 
picture [77, 78, 79]. 

Another interesting application of a photosensitive BZ system, presented 
by Rambidi and Yakovenchuk [80], is finding the shortest path in a labyrinth. 
The authors of [80] project the picture of the labyrinth onto a properly pre- 
pared Ru catalyzed system and then register the consecutive stages of chem- 
ical wave spreading throughout the image of the labyrinth’s pathway with a 
video camera. The images are processed by a computer, using a special iter- 
ative procedure which enables to establish the shortest path between given 
starting and target points. 

Another class of chemical devices which process information utilize the 
fact that a chemical pulse, which can be produced in excitable or oscillatory 
media, carries information. The area of a high concentration of a particular 
reagent may be considered as corresponding to the logical ’true” state, while 
the area of a low concentration - to the logical ” false”. Alternatively, a pulse 


present in a selected part of a system may corresponds to the logical ” 1” 
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state and if it is not there - it implies the logical ”0”. The travelling pulses 
or regular trains of such pulses are called ” chemical signals”. On such signals 
one can perform operations analogical to electronic ones. Using this idea the 
chemical reactors which work as the basic logical gates (AND, OR, NOT) 
were invented and constructed in real experiments [81, 82]. They utilize 
the BZ reaction proceeding in capillary tubes [81] or on a catalyst loaded 
membrane [82]. 

It has been found that a chemical system may perform more complex 
logical operations if it is inhomogeneous and composed of ” active” regions, 
in which reactions occur, and ” passive” areas, where some of the reagents are 
absent and so only a part of reactions proceed there. In practice it means 
that an immobilized catalyst is inhomogeneously distributed in space. The 
passive areas are those which do not contain catalyst. In such systems a 
” chemical diode”, which conducts travelling pulses only in one direction was 
constructed [83]-[86] (the structure and functioning of the chemical diode is 
presented in Section 2.4.7). The chemical diode was then used as a vital 
component of successively proposed logical gates, coincidence detectors [86] 
and memory cells, which can be written, read and erased [86, 87]. Most of 
the studies concerning such reactors present theoretical (numerical) results 
[84]-[87], but there are also a few experimental works [83, 87]. Of course, 
information processing devices of this kind are far slower than silicon ones. 
However, they may be competitive to traditional ”chips” in different fields, 


such as new kinds of 'intelligent” materials or medicaments. 
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1.4 Experimental methods. 


The numerical results presented in this thesis concern two dimensional re- 
actors, in which BZ-type reactions proceed. The considered reactors consist 
of active areas, containing the catalyst and passive areas, without catalyst. 
This section provides basic information about obtaining such systems in real 
chemical experiments. 

Typical reactions used in such experiments include the ferroin catalyzed 
BZ reaction ([19, 70, 81, 82, 83], [87]-[95]) or several types of photosensitive 
BZ reactions (catalyzed with Ce, Ru or Mn [78], [96]-[100]). 

The active areas of the system, where the reaction can proceed, are cov- 
ered with an immobilized catalyst. In case of ferroin, it may be fixed on thin 
cellulose nitrate or polysulfone membranes just by pouring the membranes 
into appropriate solutions (the ”classical” description of ferroin’s fixation 
on membranes can be found in [89]). Other cation-exchange membranes 
[70, 82, 88, 90, 91, 93, 94], [101]-[105] or mesoporous glass [83] may also be 
used. The important question is how to obtain catalyst-free pattern on an 
otherwise catalyst-covered membrane. In simplest cases the solution con- 
taining ferroin can be put on desired areas of a membrane with a fiber-pen 
[89] or printed with an ink-jet computer printer [82, 101, 102, 106], which 
provides better accuracy." The most accurate technique of fabrication of cat- 
alyst micropatterns has been recently developed by K. Suzuki, T. Yoshinobu 
and H. Iwasaki [70, 103, 104]. In their method a membrane is first covered 
with a photoresist layer on which a mask containing the desired pattern is 
printed. After exposure to light the membrane is developed and immersed in 


a ferroin solution. The ferroin molecules can then stick only to the parts of 


'Printing the catalyst with an ink-jet printer is also used in completely different sys- 


tems, eg. when desired patterns for metal deposition on a surface are needed [107]. 
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the membrane that were previously exposed to light. As the last step of the 
process the photoresist is removed by ethanol. The photolitographic tech- 
nique enables to obtain catalyst-free (passive) stripes of the width of 50m 
[70]. On the other hand, it enables to control catalyst-filled areas with 20um 
precision [103]. If a photosensitive BZ reaction is used, the desired pattern 
can be just projected onto the system by using a proper mask [78, 87]. In 
this case the reaction proceeds in the dark areas and it is inhibited in the 
illuminated ones (so they become passive). 

Travelling waves or pulses of excitation in such systems may be initiated 
very easy by touching the membrane with a silver wire [70, 82, 89, 91, 93], 
with a glass tip [105] or just by dropping the activator from a micropipet at 
the desired point of the system [108]. More sophisticated methods include the 
electrochemical method, in which a silver wire electrode present in a system 
is positively biased for a short moment (eg. 100msec), thus producing a 
chemical pulse or a wave [69, 70]. The pulses may also be initiated with an 
external electric field [109] or with a laser [94, 98]. In order to obtain a train 
of pulses with a constant frequency it is also common to use a rotating spiral, 


located in the appropriate part of the system [70, 99]. 
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1.5 The aims of this work. 


The methods of direct information processing by means of chemical reactors 
are intensively investigated, as they seem very promising. High hopes are 
associated with two dimensional reactors composed of active and passive 
areas, in which the information is carried by travelling pulses of excitation. 
Quite a large number of such reactors have already been proposed, but most 
of them come as the result of purely theoretical/numerical studies and only 
a few have been tested experimentally. Consequently, the aims of this work 


are the following: 


1. Numerical calculations oriented on finding the new geometries of reac- 
tors (ie. distribution of active and passive areas in space), which might 
be applied to direct processing of chemical signals. The structures of 
the cross junction (Section 2.3) and the simple passive barrier (Section 


2.4) have been proposed. 


2. Theoretical investigation of the phenomena which occur in proposed 
signal processing reactors. It has been found that the cross junction 
may work as a logical AND gate, a time coincidence detector of two 
pulses and a switch of the direction of propagation of a pulse (Section 
2.3). The passive barrier has the ability of systematical elimination of 
pulses from a train. Consequently, it can be used as a transformer of 


frequency of regular trains of pulses (Section 2.4). 


3. Suggesting the experiments, in which these structures and phenom- 
ena would be tested, thus proving the correctness or revealing the in- 
correctness of models used in numerical studies. Experiments on fre- 


quency transforming by a passive barrier have been made (Section 2.7). 
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Their results confirm the numerical predictions and the correctness of 
the Rovinsky - Zhabotinsky model of the ferroin catalyzed Belousov 
- Zhabotinsky reaction used here (the model is described in Section 


212) 


4. Using the simple signal processing elements the properties of which are 
already known, I can suggest more complex information processing re- 
actors, like for example a reactor that is able to add two given numbers 


together (Chapter 3). 


Chapter 2 


Results. 


2.1 The models. 


The major part of my work has been dedicated to theoretical studies of 
propagation of excitable pulses in two dimensional reactors, whose geometry 
enables processing of chemical signals in the form of single travelling pulses 
or regular trains of such pulses. The investigations have been carried out 
by numerical integration of reaction - diffusion equations describing the time 
evolution of variables appearing in selected models. Most of the results have 
been obtained for the FitzHugh - Nagumo type model (FH-N) or the Rovin- 
sky - Zhabotinsky model (R-Z). The variables of the latter model are directly 
related to concentrations of chemical spices in the ferroin catalyzed Belousov 


- Zhabotinsky reaction. Both models are presented below. 


2.1.1 The FitzHugh - Nagumo model. 


The FitzHugh-Nagumo type model used in this work [54, 55, 86] consists of 


two variables: u and v. The time evolution of their values is given by the 


27 


CHAPTER 2. RESULTS. 28 


following set of equations: 


rę = —y[ku(u — a)(u — 1) + v] + D,V7u (2.1) 
Ov 
P (2.2) 


This model is relatively easy for theoretical and numerical analysis, as it 
contains only the polynomials of the third order. It was originally developed 
to describe the propagation of pulses of excitation in nerve tissues [54, 55], 
as a simplification of the Hodgkin - Huxley model [56]. The FH-N model has 
a very good biochemical background (presented in a clear way eg. in [110]). 
The crucial disadvantage of this model is that its variables do not represent 
directly concentrations of a known reaction scheme, so the results obtained 
within this model cannot be regarded as a quantitative basis for a chemical 
experiment." Nevertheless, the model has been used as generic model for 
nonlinear chemical phenomena [57], because its dynamics reflects the basic 
features of an excitable system. 

The variables u and v in Eqs. (2.1-2.2) can be either negative or positive. 
One may say that Eqs. (2.1-2.2) describe the time evolution of a “chemical” 
system with two components. One of them (described by u) is the “activa- 
tor” (because it may speed up its own production) and the other (described 
by v) is the “inhibitor” (because it slows down the production of u). In 
this sense the values of u and v are related to the “concentrations” of the 
corresponding “reagents”. In the spatially distributed system it is assumed 
that the “inhibitor” is always immobilized (it never diffuses). The system 
(2.1-2.2) is excitable for the following parameters describing reactions: y=1, 


T=0.03, k=3.0, a=0.02 [86]. The nullclines of Eq. (2.1) without the diffu- 


1 Although there was a successful attempt to propose a microscopic reaction scheme 


corresponding to this model [58]. 
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sion term (solid line) and Eq. (2.2) (dashed vertical line identical with the v 
axis) for these values of parameters are shown in Figure 1. Solid and dashed 
arrows indicate the vector fields driving u and v respectively. The nullclines 
have only one common point u=v=0, indicating the stationary homogeneous 
solution of Eqs. (2.1-2.2). The dotted line with symbols shows a possible 
evolution of the system in the phase space after a sufficiently large perturba- 
tion. This is where the excitability of the system can be observed (note, that 
the initial perturbation does not decrease, but at first it increases).? The dif- 
fusion coefficient of the “activator” D,, equals 0.00045 because for this value 
well defined pulses propagate in the system |111, 112]. Areas of the studied 
FH-N systems where Eqs. (2.1-2.2) hold are called excitable or active areas. 

In some parts of considered structures the “inhibitor” v is absent, no 
reactions occur and only diffusion of the “activator” u is possible. These 
areas of the FH-N systems are called diffusion or passive areas. Within them 
the time evolution of u and v is described by the following set of equations 
[86]: 

Ou 


—=D),V? 2. 
Ti Vu (2.3) 


v = 0 = const. (2.4) 
where k=0.03 and D„=0.00045 [111, 112], as previously. In all the studied 
structures free diffusion of the “activator” u between active and passive areas 


is assumed. 


The response of the system described by Eqs. (2.1-2.2) without the diffusion term to 
a perturbation of v;=-0.2 is illustrated in Fig_Ola.htm on the CD-ROM (folder Figures). 
Response of the same system to two consecutive perturbations is presented in Fig_01b.htm 
and Fig _O01c.htm. In the case of Fig_01b.htm the second perturbation comes too early after 
the initial one and consequently is disappears almost instantly. In Fig_01c.htm the second 
perturbation occurs a moment later and the system is excited again. In both cases the 


second perturbation has the same amplitude of (Au=0.4, Av=0.0) 
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Figure 1. The nullclines of Eq. (2.1) without the diffusion term (solid line) 
and Eq. (2.2) (dashed vertical line identical with the v axis) for the values of 
parameters given in Section 2.1.1. Solid and dashed arrows indicate the vector 
fields driving u and v respectively. The nullclines have only one common point 
u=v=0, indicating the stationary homogeneous solution of Eqs. (2.1-2.2). The 
dotted line with symbols shows the evolution of the system in the phase space 


after a sufficiently large perturbation (vjni=-0.2). 
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One may notice that the number of parameters used in Eqs. (2.1-2.4) 
may be reduced without losing generality, which is true, as Eq. (2.1) can 
be divided by 7 and then only four independent parameters remain in the 
model. On the other hand, the form of the equations presented above is very 
convenient for manipulations, because k and a influence the threshold value 
of perturbation exciting the system, D, controls the propagation of excitable 
pulses, 7 corresponds to different time scales of the dynamics of “activator” 
and “inhibitor” and y is just a ”switch” between active and passive areas 


(y=1 for the active medium and 0 for the passive one). 


2.1.2 The Rovinsky - Zhabotinsky model. 


Another important model used in my calculations is the Rovinsky - Zhabotin- 
sky model [44, 45] of the ferroin catalyzed Belousov - Zhabotinsky reaction 
[10]. This model is based on the Field - Kórós - Noyes mechanism [42]. 
To avoid duplicating information, only the most important facts about this 
model are included in this section, because the variables, coefficients and de- 
tails about scaling can be found in Appendix A, together with the detailed 
step-by-step derivation of the reaction - diffusion equations from the chemical 
reaction scheme. 

The model consists of two variables: x and z, which correspond to dimen- 
sionless concentrations of the activator z (HBrO2) and of the oxidized form 
of catalyst z (Fe(phen)3+), where the scaling is given by Eqs. (A24-A25) 
in Appendix A. In the active regions, which contain the catalyst, the time 
evolution of the concentrations of z and z is described by Eqs. (2.5-2.6): 


dx 1 = 
5 = Tiel- 2) — (2901 87 + że (2.5) 


(2.6) 
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In the passive regions, without catalyst, the concentrations of z and z evolve 


according to Eqs. (2.7-2.8): 


ox 1.3 Tp 5 
== f f 2. 
OT al SET ve M 
z = 0 = const. (2.8) 


a, B, u and € are scaled coefficients defined by Eqs. (A28-A31) in Appendix 
A. T and p are dimensionless variables describing time and distance respec- 
tively. Their connection with physical time and distance is given by Eqs. 
(A26-A27), Appendix A. 

Eqs. (2.5-2.8) correspond to a typical experimental situation in which 
the catalyst is immobilized on a membrane, whereas the activator is in the 
solution and it can diffuse (compare eg. [19, 70, 81, 82], [87]-[91], [93, 94, 95, 
102, 113, 114)). 

In numerical calculations for the BZ system the following values of pa- 
rameters (the same as considered in [44, 45, 112], [115]-[118]) have been used: 
A=0.02 M, B=0.2 M, C=0.001 M, kj =100M~? /sec, ką=1.7-:10*M */sec, ks=107M */sec, 
kr=15M */sec, Ks=2-10 $M/sec, ky3=10 Śsec "1, q=0.5. The corresponding 
values of scaled parameters a, 8, e, p are 0.017- hę”, 0.0017- hę”, 0.1176 and 
0.00051 respectively, where ho denotes the Hammett acidity function, de- 
scribing the effective proton concentration [46, 119, 120] and it is expressed 
in mol/l. For these values of parameters the system becomes excitable if 
ho<0.9899 [118]. Here hy=0.5 has been used. 

For these values of parameters Eqs. (A26-A27) give the following con- 
stants which transform dimensionless time (7) and distance (p) into physical 
units: 


t|sec|] = 8.5- T (2.9) 


r[cm] = 2.915 - y Dx/|cm?/sec] - p (2.10) 
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where £ and r are physical time and distance. Consequently, similar transfor- 
mation constant connecting the dimensionless (v) and physical (v) velocities 


of a pulse can be calculated: 


vlem/min] = 20.58 -  Dx/|cm2/sec] -v (2.11) 


In the last two equations Dy is the diffusion constant of the activator X. How- 
ever, this constant strongly depends on the medium the reactions proceed in. 
In the aqueous solution it is of the order of 10~°cm?/sec [45, 46, 84, 118, 119] 
whereas for a reaction in a gel it may be reduced by two orders of magnitude 
[61, 84]. In order to make the results more general, all distances and veloc- 
ities concerning the R-Z model are presented in double form: dimensionless 
and as the function of the ratio of diffusion constants Da where the value 
of Dx, corresponds to a particular choice of the diffusion constant: Dx, =1: 
105cm? /sec [45, 118]. The second number allows to see more clearly the 
real spatial and temporal scale of the considered processes. 

The nullclines of Eq. (2.5) without the diffusion term (solid line) and 
Eq. (2.6) (dashed line) for these values of parameters are shown in Figure 
2. The dotted line with symbols shows a possible evolution of the system in 


the phase space after a sufficiently large perturbation.® 


The coordinates of intersection of nullclines in Fig. 2 give the stationary 


3Like in the case of the FH-N model, there are additional animated figures illustrating 
nullclines of the system described by Eqs. (2.5-2.6) without the diffusion term (CD-ROM, 
folder Figures). The response of such a system to a perturbation of x;=0.1 is illustrated 
in Fig_02a.htm Response of the same system two consecutive perturbations is presented 
in Fig_02b.htm and Fig_02c.htm. In the case of Fig02b.htm the second perturbation 
comes too early after the initial one and consequently is disappears almost instantly. 
In Fig_02c.htm the second perturbation occurs a moment later and the system is excited 


again. In both cases the second perturbation has the same amplitude of (Ar=0.1, Az=0.0) 
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Figure 2. The nullclines of Eq. (2.5) without the diffusion term (solid line) and 
Eq. (2.6) (dashed line) for the values of parameters given in Section 2.1.2. Solid 
and dashed arrows indicate the vector fields driving z and z respectively. The 
nullclines have only one common point (%sa=7.283-10-4, %sa=1.060-10-2) indicat- 
ing the stationary homogeneous solution of Eqs. (2.5-2.6). The dotted line with 
symbols shows the evolution of the system in the phase space after a sufficiently 


large perturbation (2j,;=0.1). 
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values of z and z in the active (excitable) area. They read: 
Esa = 7.283 - 1074 (2.12) 


Zsa = 1.060 - 107? (2.13) 


(which is the stationary solution of Eqs. (2.5-2.6)). The stationary concen- 

trations in the passive area (the stationary solution of Eqs. (2.7-2.8)) are 
given by: 

Tsp = 5.099 - 1074 (2.14) 

ży =0 (2.15) 


Free diffusion of the activator x is assumed between active and passive areas. 
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2.2 Propagation of single pulses of excitation 
through a passive gap. 


I. Motoike and K. Yoshikawa considered two active areas separated by a stripe 
of the passive one [86]. A pulse of excitation propagating in one of these 
areas may excite the other active area if the passive stripe is narrow. The 
maximum width of the passive stripe for which such excitation still occurs 
is called the penetration depth. They found [86] that the penetration depth 
depends on the geometry of the junction and on the direction of propagation 
of incident pulses and it is maximal for plain pulses travelling in the direction 
perpendicular to the barrier. In [86] Motoike and Yoshikawa describe their 
results obtained with an Euler explicit method for the FitzHugh - Nagumo 
model, the same and with the same values of parameters as considered here. 
According to [86], if the penetration depth for a single pulse with the wave 
vector perpendicular to the barrier is denoted by d..ry_n, then for a pulse 
travelling parallel to the barrier it is only 0.94-d. ry_N. This result has been 
confirmed for the purpose of this study with an implicit method based on the 
Crank - Nicolson discretization of the Laplace operator [121].* It has been 
found [86, 116] that d.ry—n0.163. 

Similar calculations have been carried out for the Rovinsky - Zhabotinsky 
model [115, 116]. The penetration depth for a single pulse travelling in the 
direction perpendicular to the barrier is de p_z¥3.295 (0.03037: DCM), 
while for a pulse travelling parallel to the barrier it is about 3.287 (0.03029- Dzzem), 
which is 0.995-de r-z. Note, that this result implies, that the penetration 


4This method has been used to obtain most of the results presented in this work; if 
other algorithms have been applied, it is indicated in the text. The formulas corresponding 
to implicit integration schemes mentioned here as well as the explicit Euler method can 


be found in Appendix B. 
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depth for “perpendicular” and “parallel” pulses differ by 0.5% only, while 
for the FitzHugh - Nagumo model it is 6%. Creating a stripe of passive 
field with higher than 0.5% accuracy in its width is a hard, but achievable 
task from the experimental point of view. One may apply for example the 


photolithographic methods for the introduction of catalyst [70, 103, 104].° 


5 At first the BZ system was studied for ho=0.97. The penetration depths for pulses 
travelling perpendicularly or parallel to the passive stripe for this value of hg are almost 
the same, as they differ by less than 0.05%. Therefore it seems to be very difficult to make 
a passive stripe impenetrable for a ” parallel pulse” and penetrable for a ” perpendicular” 
one in the laboratory. Decrement of ho stabilizes the stationary solution of Eqs. (2.5-2.6) 
[115]. This results in a faster relaxation of the active medium and remarkably increases 
the asymmetry in the penetration depth. From the chemical point of view smaller ho 
corresponds to lower concentration of H* in the system. The results given in text have 


been obtained for ho=0.5. 
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2.3 The cross junction and its properties. 


Having in mind the results of the previous section (Section 2.2) concerning the 

FH-N model one may notice that a layer of the width dpy_y (0.94-d. ry_n <dry_n <de PH-N) 
is transparent for a pulse propagating perpendicularly, but impenetrable for 

a pulse propagating parallel to it. This observation opens the way to con- 

struct a simple coincidence detector of chemical signals. 

Let us consider a plain (two dimensional) system with geometry shown 
in Figure 3, where the brighter areas correspond to the active (excitable) 
medium and the dark stripes - to the passive (diffusion) one. The excitable 
and diffusion areas in the device shown in Fig. 3 form two perpendicular 
“channels”, in which travelling pulses can propagate. To simplify the de- 
scription that follows let us assume that the excitable areas of our cross 
junction are numbered as in Fig. 3. 

If the width of the stripes of the diffusion field dpy_y is between 0.94- de FH-N 
and d, then the cross junction presented in Fig. 3 may have interesting prop- 


erties. It can be expected that: 


i A single pulse propagating in one of the channels is able to propagate 
along its own channel, but it cannot move “sideways” (because it cannot 
cross the stripe of the diffusion field parallel to the direction of its 


motion, as dpH-N>0.94:d. FH-N); 


ŚThe asymmetry of the penetration depth within the FH-N model is nicely illustrated 
in Fig 03a.htm and Fig_03b.htm on the CD-ROM (folder Figures). Brighter colors cor- 
respond to higher concentration of the “activator” u. The passive stripe is blue and the 
active areas are red. The same color coding is used for all the other animations mentioned 
later in text. In Fig_03a.htm a plain pulse of excitation gets through a barrier perpendic- 
ular to the direction of its motion, while the same barrier cannot be crossed by a pulse 


travelling parallel to it (Fig_03b.htm). 
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Figure 3. Geometry of the cross junction: brighter areas correspond to the 
excitable field and the dark stripes - to the diffusion one. The excitable areas 


are numbered to simplify the description in text. 


http://rcin.org.pl 
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ii The pulse is able to propagate through the junction of two channels 
(because it is able to cross the stripe of the diffusion field perpendicular 


to the direction of its motion, as dpry_N<d-FH-N). 


Of course, properties (i) and (ii) are formulated for the FitzHugh - Nagumo 
model, but analogical construction and explanations are valid for the Rovin- 
sky - Zhabotinsky model.” 

If there are two pulses of excitation propagating in perpendicular channels 
of the cross junction presented in Fig. 3, they may coincide and interact in the 
common central part of the device (area 5 in Fig. 3). This interaction may 
be different, depending on the time difference between the incoming signals. 
Thus, the detection of coincidence of signals from individual channels should 
occur in the central area of the cross junction. It is interesting and useful 
to obtain the full characteristics of the behavior of the cross junction with 


respect to the time difference between the incoming pulses. 


2.3.1 Results for the FH-N model. 


The properties of the cross junction within the FH-N model have been thor- 
oughly studied for passive stripes of the width dpy_Nn=0.16 [111]. The pulses 
have been initiated on the border by decreasing the value of the inhibitor v 
inside the channel to v;=-0.2. The answer of a homogeneous system with 
respect to such excitation is shown in Fig. 1 (Section 2.1.1). The stationary 
shape of a pulse is obtained after t.prH_Nn 2.50 and the velocity of free 
propagation of the pulse is vs py—_n 0.77. Thus after covering the distance 


of 2 units the pulse retains its stationary form. 


Properties (i) and (ii) for the FH-N model are presented in Fig_03c.htm on the CD- 
ROM (folder Figures). All passive stripes in Fig_03c.htm have the same width of about 
0.96- de, FH-N 
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Figure 4 presents a typical stationary shape of a travelling pulse in terms 
of the variable u. Brighter areas correspond to the higher values of u. Thin 
stripes of the diffusion field are located between the pairs of dashed lines 
(the stripes are not filled black, so that the picture is more legible). One can 
see characteristic “whiskers” of the pulse, slightly sticking outside the stripes 
of the diffusion field. The collision of these “whiskers” explains part of the 
characteristics of the cross junction given below (points A and C). Such sta- 
tionary shape of a pulse does not depend on the mechanism of initialization. 

If the width of the excitable field inside the channel (de FH_N) is large 
enough (de rH-N>0.40) we do not observe its influence on the amplitude and 
velocity of a pulse. In very narrow channels (de ru-n<0.10) both amplitude 
and velocity are remarkably smaller than those characterizing a pulse in a 
wide channel. In particular, in a channel characterized by de ru_N=0.10 the 
amplitude of the pulse of u (” activator”) is by 2% smaller than for a pulse 
propagating in a wide channel (uma1=0.93 has been observed instead of 0.95 
seen in wide channels). The velocity of such a pulse drops by 13%, from 0.77 
to 0.67. In the narrowest channels studied (de ru-N=0.02) the amplitude of 
the pulses of u iS umaz =0.75 and their velocity is only 0.39, but the pulses 
still propagate. 

The features of a pulse mentioned above are important for its propagation 
through the junction. The data concerning penetration of diffusion filed 
given by I. Motoike and K. Yoshikawa in [86] applies to pulses propagating 
in wide channels. If the width of the channel is smaller than 0.70, then the 
penetration depth for transverse propagation decreases. The width of the 
excitable field inside each channel is also important for a pulse to propagate 
through the junction, because the pulse needs enough space between the two 


stripes of the diffusion field on its way to cross both of them. In other words: 
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Ba 1.0 2.0 30 40 50 60 70 


Figure 4. Typical stationary shape of a travelling pulse in terms of the 
variable u. Brighter colors correspond to higher values of u. Thin stripes of 
the diffusion field are located between the pairs of dashed lines (the stripes 
are not filled black, so that the picture is more legible). 


http://rcin.org.pl 
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to obtain the propagation of a pulse from the area 4 to the area 6 in Fig. 
3 - the area 5 cannot be too small. It has been found [111] that in order 
to obtain the junction with features (i) and (ii) the inner diameter of the 
channel de run should not be lower than ca. 0.80. 

All the results presented below are obtained by integrating Eqs. (2.1-2.4) 
numerically. The calculations have been carried out with the implicit method 
[121]. The values of parameters given in Section 2.1.1 have been used, as well 
as de ru_N=l.02. The time step of the integration dt=0.005. Three sizes 
of the square the junction is built on have been considered (5x5, 7x7 and 
8x8), covered with square lattices of 250x250, 350x350 or 400x400 points 
respectively. Thus, the space step in all those calculations dr=dy=0.02 re- 
mains constant. The results for different system’s sizes are consistent. No 
flux boundary conditions have been used at the border of the square. 

To find the time characteristics of the cross junction the following tech- 
nique has been applied. In each of the input channels (areas 4 and 8) one 
input signal (a pulse) is initiated. The pulses are produced on the border of 
the square. For various experiments the channels are located in different parts 
of the square, so that the pulses travel over various distances before reaching 
the cross junction. This way, due to their constant velocity (vs rH-N 0.77), 
one is able to adjust the time difference At between their arrival at the junc- 
tion. Depending on At seven types of the junction’s response have been 


obtained. They are presented in Figures 5-11. The figures marked as a)-f) 


SNote, that At can be specified only with finite resolution of 
da /vs FH-N=0.02/0.77=0.026, as the distance between the place of initiation of a 
pulse and the cross junction is controlled up to the space step of the grid used (0.02). 


This makes the ”spectrum” of At discontineous (see points A.-G. below). 
Each of Figs. 5-11 has its color, animated version under the name like Fig_05a.htm, 


Fig_06a.htm, ... Fig_lla.htm (CD-ROM, folder Figures). 
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show the time evolution of the activator in the junction (the exact timing of 
snapshots a)-f) is given in figure captions). 

Due to the symmetry of the system it is enough to restrict the analysis 
to the case in which the first pulse is initiated in the area 4 and the second 
- in the area 8. The behavior of the cross junction may be classified in the 


following way: 


A. For a very small time difference between incoming signals (At<0.26) the 
cross junction acts as the “AND” logical gate. The output signal is 
produced in the area 7 between the two input channels (4 and 8) if and 
only if there are two input signals arriving at the junction within the 
given time difference. The output signal probably appears when the 


“whiskers” of the two travelling pulses meet (cf. Fig. 5). 


B. For the time difference in the range 0.28< At<3.45 only the earlier signal 
survives and propagates within its own channel - area 6. Here the cross 
junction works (as expected) as a coincidence detector in which the 
earlier pulse switches the central field into the refractory regime and 


blocks the propagation of the second signal (cf. Fig. 6). 


C. When 3.48< At<3.79 the earlier signal survives in its own channel (area 
6) and in addition a new output signal is produced in the lower right 
hand side corner of the device (area 9). The new output signal appears 
due to the same mechanism as in the case A, ie. due to the coincidence 


of the “whiskers” (cf. Fig. 7). 


D. For 3.81< At<4.22 the evolution is the same as in the case B: only the 
earlier incoming signal survives and propagates in the area 6 (cf. Fig. 


8). 
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Figure 5. The “AND” logical gate for At=0.00. The output signal is 


produced in the area 7 between the two input channels (4 and 8). The 


consecutive snapshots correspond to a) t=0.5, b) t=2.5, c) t=3.5, d) t=4.0, 
e) t=4.5, f) t=5.0. Such behavior is observed for Aże[0.00, 0.26]. 
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Figure 6. Only the earlier signal survives (within its own channel - area 6) 
for At=1.29. The consecutive snapshots correspond to a) t=0.5, b) t=2.5, c) 
t=4.0, d) t=4.5, e) t=5.5, f) t=7.5. Such behavior is observed for AtE|0.28, 
3.45]. 
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e) 


0. 0. x 
bo 1.0 2.0 3.0 4.0 5.0 6.0 bo 1.0 2.0 3.0 4.0 5.0 6.0 


Figure 7. The earlier signal survives in its own channel (area 6) and in 
addition a new output signal is produced in the area 9 for At=3.61. The 
consecutive snapshots correspond to a) t=1.0, b) t=4.0, c) t=6.0, d) t=7.0, 
e) t=8.0, f) t=10.0. Such behavior is observed for Ate[3.48, 3.79]. 
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Figure 8. For At=3.87 only the earlier incoming signal survives (area 6). 
The consecutive snapshots correspond to a) t=1.0, b) t=4.0, c) t=6.0, d) 
t=7.0 e) t=8.0 f) t=9.0. Such behavior is observed for Ate[3.81, 4.22]. 
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Figure 9. For At=4.25 the earlier incoming signal follows the way within its 
own channel (area 6) while the other signal is switched from the original path 
and follows the first signal (area 6). The consecutive snapshots correspond 
to a) t=4.0, b) t=6.0, c) t=8.0, d) t=8.5, e) t=10.0, f) t=11.0. Slides c) and 
d) present ” zoomed” images of the central field of the junction, in order to 


show in details the most interesting part of the evolution in this area. 
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d) 


e) 


F 0. 
00.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 bo 1.0 2.0 3.0 4.0 5.0 6.0 7.0 


Figure 10. For At=4.64 the earlier signal goes unchanged through the 
junction, while the other signal initiates output pulses both within its own 
channel (area 2) and in the horizontal channel (area 6). The consecutive 
snapshots correspond to a) t=4.0, b) t=7.0, c) t=8.5, d) t=9.0, e) t=10.0, 
f) t=11.5. Such behavior is observed for Ate[4.28, 5.41]. 
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e) 


0. 0. a 
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Figure 11. For At=5.67 the pulses do not “feel” each other any more. 
Each of them simply propagates in the channel it was initialized in (areas 6 
and 2 respectively). The consecutive snapshots correspond to a) t=4.0, b) 
t=6.0, c) t=7.0, d) t=9.5, e) t=10.0, f) t=12.0. Such behavior is observed 
for At>5.43. 
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E. For Atx4.25 the earlier incoming signal propagates within its own chan- 
nel (area 6), while the other signal is switched from the original path 
and follows the first signal (cf. Fig. 9). This behavior is caused by the 
fact that the second pulse finds the central area of the junction relaxed 
only on the left hand side. Consequently, further propagation of the 
pulse within the junction starts there and the geometry of the reactor 
makes the pulse switch its original direction. It can be seen in Figs. 9c) 


and 9d), which present the ”zoom” of the central filed of the junction. 


F. For 4.28< At<5.41 the earlier signal passes unchanged through the junc- 
tion, while the later signal initiates output pulses both within its own 
channel (area 2) and in the horizontal channel (area 6). The explana- 
tion of this case is similar to the one given it the previous point. The 
only difference is that here the second pulse enters the central area of 
the junction a bit closer to the right than in the case E, due to a better 


relaxation of the active medium in the central area (cf. Fig. 10). 


G. For A£>5.43 the pulses do not “feel” each other any more. Each of them 
simply propagates in the channel it was initialized in (areas 6 and 2 
respectively). This is because the central area of the junction is fully 


relaxed when the second pulse appears (cf. Fig. 11). 


2.3.2 Results for the R-Z model. 


Interesting results obtained within the FitzHugh-Nagumo model motivated 
me to study the propagation of excitable pulses in the cross junction geome- 
try, but using a model which can be directly linked with chemical reactions, 
in order to give grounds to verify the results experimentally. For this purpose 


the Rovinsky - Zhabotinsky model has been applied. 
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The system of partial differential equations (2.5-2.8) has been solved 
numerically using the explicit Euler technique with a constant time step 
dr=0.001. The pulses have been initialized at the boundaries of the square, 
inside the channels (areas 4 and 8 in Fig. 3) by increasing the concentration 
of z to 0.1. In the majority of calculations the square has been covered with 
a grid of 320x320 points. At the beginning the concentrations of z and z 
in all active and passive areas correspond to their stationary states Zsa, Zsa 
and Zsp, Zsp respectively (Eqs. (2.12-2.15)). On the lower and left hand side 
borders of the system no flux boundary conditions have been used. On the 
other hand the values of z and z have been fixed as equal to the correspond- 
ing stationary states (Zsa, Zsa aNd Tsp, Zsp respectively) at the upper and right 
hand side borders. This way the pulses arriving there may freely ”leak out” 
of the system (they do not accumulate in the neighborhood of the borders) 
and any influence of the system’s boundaries on the time evolution can be 
avoided. 

Passive stripes which are dp_z=3.291 (0.03034- pcm) wide have been 
considered (this value is between 0.995:d.p_z and de-z given in Section 
2.2). From the calculations the minimum width of a signal channel has been 
found. For the selected width of the passive stripes the pulse initiated (for 
example) in the area 8 always arrives to the central area 5. However, if 
the central field is too small then the amplitude of the activator is not high 
enough to excite the opposite area 2. It has been found that the minimum 
inner width of the signal channel for which pulses can propagate through the 
junction is about 10.526 (0.097029. pcm). 

The typical velocity of the pulse in our system is about 4.128 (0.2686- Dzzem/min), 
which is a reasonable value for a chemical wave in the BZ reaction [46, 102, 


Fis; M0, 122 122]. 
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Knowing the width of a semitransparent passive barrier and the width of 
an active channel one can build a cross junction of the excitable BZ system 
which should work similarly to that for the FitzHugh - Nagumo dynamics. 
The pulses initiated in the horizontal channel (in the area 4 in Fig. 3) should 
be able to get through the junction (area 5 in Fig. 3) and arrive at the other 
side of the device (area 6 in Fig. 3) without escaping from the channel (to 
the areas 1-3 and 7-9 in Fig. 3) and due to the symmetry of the device the 
pulses in the vertical channel should behave in the same way. What is the 
interaction between pulses coming from perpendicular directions? At the be- 
ginning the system is in its stationary state. In calculations two sizes of the 
square on which the junction is built have been considered: 198.13x198.13 


and 263.3x263.3 (corresponding to 1.8264- 5ź cm and 2.427- Džem), for 
0 (0) 


which square lattices of 240x240 or 320x320 points respectively have been 
used. Thus, the space step in all those calculations dp=0.82 remains con- 
stant. Two channels of the width of 66.648 (0.61436- Dz m) or 132.473 
(1.22114- Dzem) for smaller and larger systems respectively are placed 
symmetrically on the square. The results obtained for different grids are 
consistent. The first pulse is initiated on the left hand side border of the 
area 4 in Fig. 3. After a time Ar another pulse is initiated on the lower 
border of the area 8. Because the distance from both borders to the junction 
is the same, the time difference between the arrival of the first and the second 
pulse at the junction is also Ar. 

It has been observed that if the second pulse arrives at the central field 
earlier than ATmin=107 (910sec) after the first one, it is stopped at the central 
field. On the other hand, if AT > Armaa=115 (978sec) both pulses follow 
their original signal paths without interaction. 


Although it is apparently possible to build a coincidence detector based 
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on the cross junction of excitable BZ systems with the dynamics given by the 
Rovinsky - Zhabotinsky model, the other interesting features of the junction, 
described in Section 2.3.1, do not occur in the R-Z model. In particular, the 
” AND” logical gate behavior, when two pulses coincide, producing an output 
pulse between the channels (in the areas 7 or 9 in Fig. 3) has never been 
observed. Neither the ” chemical switch”, when the second pulse changes the 
direction of propagation and follows the first one [111]. 

In order to see more clearly why the ” chemical switch” does not work for 
the BZ system, calculations for the asymmetric junction presented in Figure 
12 have been performed.'® Here a square of 263.3x263.3 (2.427: Džem), 
covered with the grid of 320x320 points has been used. The first pulse is 
initiated in the horizontal channel (area 4 in Fig. 3). The channel is located 
at the distance 56.363 (0.51956- Dzem) from the lower boundary; its width 
is 99.560 (0.91775- Dzem). After the time Ar=109.0 (927sec) another 
pulse is initiated in the vertical channel (area 8 in Fig. 3). This channel 
is located at the distance 56.363 (0.51956- DZ ccm) from the left hand side 
boundary; its width is 152.220 (1.40317- Bx om). One can see that when 
the second pulse gets to the central area of the junction, it propagates in both 
vertical and horizontal directions. The propagation to the right fails, while 
the propagation upwards continues and, what is very important, the pulse is 
able to pass from its own channel (area 2 in Fig. 3) to the upper right hand 
side part of the system (area 3 in Fig. 3). Comparison of Figs. 12d) and 12e) 
shows that the pulses trying to get from the area 5 to the area 6 and from 
the area 2 to the area 3 ” attack” the stripe of the diffusion field both at the 


same angle, so the geometrical factor is not the reason for the propagation 


10 Animated version of this Figure (Fig_12a.htm) can be found on the CD-ROM (folder 
Figures). 
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Figure 12. Two pulses in an asymmetric cross junction of systems with the 
Rovinsky - Zhabotinsky dynamics. Figures show the concentration of activator for 
a few selected moments in the case the second pulse of excitation arrives Ar=109 
after the first one. The brighter areas correspond to higher values of z. a) T=54 
- the first pulse passes through the horizontal channel; b) r=114 - the second 
pulse approaches in the vertical channel; c) r=134 - the second pulse propagates 
asymmetrically in the central area; d) T=150 - although the geometrical condition 
is fulfilled, the pulse cannot pass into the area 6, but it gets into the area 2; e) 


rcin.O 
pulse gets into the area 3. Such behavior is observed for Are(107, 115). 


T=170 - asymmetric propagation Ried: seecond ese in the area 2; f) r=184 - the 
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failure. The failure occurs apparently because the active medium in the area 
6 is still not relaxed behind the first pulse. A number of calculations for a 
set of smaller time shifts Ar has been performed (until the second signal is 
not able to enter the junction at all for AT=ATmin=107 (910sec)) as well as 
for the larger ones (until the second signal propagates through the junction 
and does not ” feel” the influence of the first signal at all for AT=ATmaz=115 
(978sec)) and the second signal turning right from the area 5 to the area 6 
has never been observed. Consequently, the ”chemical switch” behavior in 
the Rovinsky - Zhabotinsky model has not been found, despite the favorable 
geometry of the system. 

The difference with respect to the FitzHugh - Nagumo model may be 
explained as follows: in the Rovinsky - Zhabotinsky model the distance at 
which the catalyst (z) relaxes is by order of magnitude larger than the one 
for the activator (x), whereas in the FitzHugh - Nagumo both u and v relax 
at the same spatial scale. 

Looking at Figs. 12e) and 12f) one can see that the simplest idea how to 
make the second pulse turn right is to direct the pulse excited in the upper 
right hand side corner of the system (area 3 in Fig. 3) into an additional 
channel. 

This has been done in the computer experiment illustrated in Figure 13.1! 
The size of the square, grid size, width of the diffusion stripes and the time 
shift between the signals are the same as in Fig. 12. The main channels 
of identical width of 132.473 (1.22114- Dzem) are placed symmetrically 
in the middle of the square (at the distance 66.236 (0.61057-,/2%cm) from 


DX 


each border). An extra channel of the diameter 33.735 (0.31097-, /2* em) is 


Dx 


11 Animated version of this Figure (Fig_13a.htm) can be found on the CD-ROM (folder 
Figures). 
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Figure 13. Two pulses in a symmetrical cross junction of systems with the 
Rovinsky - Zhabotinsky dynamics with an additional channel. Figures show con- 
centration of activator for a few selected moments in the case the second pulse of 
excitation arrives Ar=109 after the first one. The brighter areas correspond to 
higher values of x. a) T=54 - the first pulse passes through the horizontal channel; 
b) r=114 - the second pulse approaches in the vertical channel; c) r=140 - the 
second pulse propagates asymmetrically in the central area; d) r=160 - the pulse 
is stopped on the right hand side boundary of the vertical channel; e) T=170 - 
asymmetric propagation of the SPA pule in t $I area 2; f) Tr=182 - the pulse 


cin.org.p 
gets into the area 3 (as in Fig. 12f), but also generates a pulse in the extra channel. 
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located in the upper right hand side corner of the system. In Fig. 13f) one 
can see that the pulse is in fact excited in that extra channel. The results 
obtained suggest that this type of behavior can be observed (at least) for 
AT € [108, 111] [918sec, 944sec]. 

Another way of building a ” chemical switch” is to speed up the relaxation 
of the active chemical medium in the area 6 (Fig. 3) behind the first pulse. 
One can do it easily by increasing the diffusion coefficient in the area 6. 
Now pulses propagate here faster and the refractory time of the medium is 
shorter. The idea is demonstrated by results of another calculations presented 
in Figure 14.1? The geometry is the same as in Fig. 13, except there is no the 
extra channel. Instead, the diffusion coefficient of the activator is doubled 
in the area 6, comparing to the rest of the system. The first pulse travels in 
the horizontal channel (areas 4-6 in Fig. 3) and speeds up remarkably after 
crossing the junction (in the area 6 in Fig. 3). After Ar=109.0 (927sec) a 
pulse in the vertical channel is initiated. This second pulse gets through the 
junction (within its own channel - area 2 in Fig. 3), but it also follows the 
first pulse, as desired. As in the previous calculations (Fig. 13) this behavior 
can be observed (at least) for Ar € [108, 111] [918sec, 944sec]. 

The calculations have shown that the width of signal channels does not 
have significant influence on the behavior of the cross junction of the Rovin- 
sky - Zhabotinsky systems (if only the channels are wider than the minimum 
0.9703mm). The system for which the signal paths are three times wider 
than those shown in Figs. 13 and 14 has been investigated (thus, the signal 


channels have been 3.66cm wide). It has been found that for this system the 


12 Animated version of this figure (Fig_14a.htm) can be found on the CD-ROM (folder 
Figures). In addition, Fig_14b.htm proves that in the last case a pulse propagating in the 
vertical channel does not ”leak out” of it, despite the doubled diffusion coefficient in the 


area 6. 
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Figure 14. Two pulses in an asymmetric cross junction of systems with the 
Rovinsky - Zhabotinsky dynamics, in the case of faster diffusion in the area 6. 
Figures show concentration of activator for a few selected moments in the case 
the second pulse of excitation arrives AT=109 after the first one. The brighter 
areas correspond to higher values of x. a) T=54 - the first pulse passes through 
the horizontal channel; b) T=114 - the second pulse approaches in the vertical 
channel; c) T=146 - the second pulse propagates asymmetrically in the central 
area; d) T=158 - the pulse is doubled: it turns right into the area 6 and follows its 


cin.O 
the areas 2 and 6; f) T=176 - the pulse from the area 6 excites a pulse in area 3. 


original path (area 2); e) T=164 nież rę Pacppeation of the second pulse in 
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pulses which arrive earlier than Ar=100 after the first one are stopped. Such 
result might have been expected as the ”transparency” of the central field is 
related to the time the system spends in the refractory regime, which hardly 
changes with size. On the other hand, the interval of times for which the 
system works as a signal switch becomes slightly longer for larger systems. 
This is related to the direction of propagation of the second pulse in the 
central area. It has been observed that for three times wider signal channels 
(3.66cm) the junction works as shown in Fig. 12, even if the second pulse 


arrives AT=120 after the first one. 


2.3.3 Discussion. 


It has been demonstrated that a cross junction of active and passive fields 
reveals interesting properties. These properties include coincidence detection 
(the pulses coming within times shorter than a characteristic one are stopped) 
and switching the direction of propagation of a chemical pulse. The features 
of the cross junction have been studied for the FitzHugh - Nagumo and the 
Rovinsky - Zhabotinsky dynamics. They can be realized in a junction of 
excitable systems of Belousov - Zhabotinsky reaction. 

It has been found that in order to construct a working cross junction 
based on the Belousov - Zhabotinsky system it is necessary to create passive 
stripes with a very high accuracy. If the diffusion constant of the activa- 
tor corresponds to the value characteristic for an aqueous system (Dx, =I: 
10-°cm?/sec) then the stripe without catalyst should be 0.3034mm wide 
with only 1% of tolerance. Such a high precision can be achieved with newly 
developed photolithographic technique of introducing catalyst on the surface 
[70, 103, 104]. 


The most interesting feature of the cross junction is its ability to switch 
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the direction of an incoming signal. The effect has been discovered for a 
symmetrical cross junction (the signal channels are of equal width) and the 
FitzHugh - Nagumo dynamics. Similar effect for the BZ system has not 
been observed, which seems to be connected to its slower relaxation. It 
has been found that a switch of a chemical signal can be constructed if the 
cross junction of excitable BZ systems is not symmetrical. Two working 
schemes, in which the asymmetry is related to the diffusion coefficient or to 
the existence of another output channel have been studied. 

It is worth noticing that in the case of Belousov - Zhabotinsky reaction 
such properties as the penetration depth strongly depend on the parameters 
of the model. Therefore the comparison between calculations and carefully 
prepared experiments for a junction of passive and active areas may be a 


good verification of the model and the values of parameters used. 
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2.4 Propagation of regular trains of pulses 
through a passive barrier. 


This section is concerned with a simple two dimensional reactor consisting 
of an infinite plane of the active (excitable) medium with a stripe of the 
passive medium (a passive barrier). Unlike in the previous section, instead 
of single pulses of excitation - regular trains of such pulses arriving at the 
barrier with a certain frequency are now considered. Such a train of pulses 
is called a chemical signal. Pulses in the form of planar waves, travelling 
perpendicularly to the stripe are considered. Such a problem is symmetrical 
in the direction perpendicular to the stripe, so it can be described as one 
dimensional. Therefore in the calculations the system’s evolution has been 
modelled on a (one dimensional) interval. 

A scheme of the investigated system is shown in Figure 15. In the bottom 
part of Fig. 15 the black lines correspond to the active areas, which are 
separated by a passive barrier (the gray line). There are free flow boundary 
conditions between passive and active media and no flux boundary conditions 
at both ends of the interval. Initially both active and passive areas are in 
their stationary states. Pulses of excitation are initiated at the left end of 
the interval and they travel to the right, coming across the passive barrier 
on their way. Initiation of pulses is performed regularly at times ktp, for 
selected tp >0 and k=1, 2, 3,...,kmax- In the following the "input signal” (a 
train of pulses arriving at the barrier) and the ”output signal” (a train of 
pulses going away after crossing the barrier) are distinguished. The input 
signal frequency is defined as f,=1/t,. If the passive barrier is narrow, it is 
transparent to all the incoming pulses; if it is wide - it is impenetrable for 


any of them. However, it turns out [116, 117, 124] that for a certain range 
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Figure 15. A scheme of the system studied. It is represented by one di- 


mensional interval, shown in the bottom part of the figure; the active areas 


are black and the passive barrier is gray. The barrier is located between grid 


points nı and ng. Concentrations of reagents are observed at grid points 4 


(before the barrier), iz (just behind the barrier) and iz (far behind the bar- 


rier). The upper part shows a snapshot of a train of pulses propagating in 


the system (incoming signal, on the left hand side and outgoing signal, on 


the right hand side). Note that the signals have different frequencies. 
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of f, and barrier’s width d the passive gap transforms the frequency of the 
input signal. It means that the frequency of the output signal fe (observed 
behind the barrier) is a fraction of f,, because some of the incoming pulses 
are stopped at the barrier, while the others get through (f,/f,€[0, 1]). This 
fact is schematically shown in the upper part of Fig. 15. 


2.4.1 The numerical technique. 


The time evolution of the signal is studied by numerical integration of the 
reaction-diffusion equations in the active areas and within the barrier (Eqs. 
(2.1 - 2.4) for the FH-N model or (2.5 - 2.8) for the R-Z model). The system, 
represented by an interval of the length /, is divided into n parts by n+1 
points of a grid, including both ends (cf. the bottom part of Fig. 15). The 
distances between the grid points dl may be different; by selecting a fine 
grid around the barrier it is possible to increase the accuracy of calculations, 
while a crude grid far away from the barrier saves computer time without 
affecting the accuracy. The barrier is located between grid points m, and m 
(<n, <n<«n, cf. Fig. 15). It means that the evolution of the system at all 
grid points 7 €[0,m1]U|n2,n] is given by the set of reaction-diffusion equations 
corresponding to the active medium (Eqs. (2.1 - 2.2) for FH-N model or 
Eqs. (2.5 - 2.6) for FH-N model), and the equations describing the passive 
medium (Eqs. (2.3 - 2.4) or (2.7 - 2.8) respectively) give the evolution at all 


grid points 7€(n,,n2). The barrier’s width is estimated as: 
dz (na — ni — 1) $ dl (2.16) 


where dh is the distance between neighboring grid points (in the calculations 
dh is constant in the barrier and around it). 


The concentrations of reagents of interest have been calculated using the 
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implicit method already mentioned with the distance between neighboring 
grid points (dl) as the space step of numerical integration. In the compu- 
tations two types of the space grids have been applied . In case of equally 
spaced (uniform) grids, used in part of calculations, the whole interval / has 
been divided into n equal parts (in this case dl=l/n). 

In the other calculations adaptive grids have been used, for which the 
passive barrier and its neighborhood (so the part of the system which has 
the most influence on signal’s transformation) has been covered with a fine 
grid and at both ends of the interval a longer distance between grid points 
has been used.'? Two kinds of adaptive grids have been applied. One of 
such grids is formed by n=800 grid points and 30 of them are placed within 
the barrier (the barrier is located between n,=400 and ng=431). The time 
evolution in the active areas surrounding the barrier is also calculated using 
the same fine grid as for the barrier di=dl/8. 40 grid points on each side 
of the barrier, with the distance dl, between them are considered. Next 
there are 40 grid points of each side with the distance 2-dl, and yet another 
40 points with the distance 4-dł4. The reaction-diffusion equations in the 
remaining part of the system are solved with a crude grid of dl (cf. Figure 
16). In order to test the numerical stability, another adaptive grid has also 
been used, for which the total number of grid points is n=870 and in this case 
the passive barrier, located between n,=420 and mm=481, is covered with 60 
points of dl,=dl/16. In the active areas located on both sides of the barrier 
the grid is fine (dl/16, 40 points) and next it increases as follows: dl/8 - 
20 points, dl/4 - 40 points and dl/2 - 40 points. The rest of the system is 


covered with the space step dl. 


13The formulas corresponding to the implicit method of numerical integration on adap- 


tive grids are given in Appendix B. 
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Figure 16. A scheme of the adaptive grid for n=800. Black intervals corre- 
spond to active areas, the gray one marks the barrier. The distance between 
grid points is given above the line, the number of grid points separated by a 


particular dl is given below the line. 
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The values of activator and inhibitor are recorded at indicators located 
at the grid points 4 (before the barrier), % (just behind the barrier) and is, 
located far behind the barrier, as shown in the bottom part of Fig. 15. By 
comparing the time evolutions at 4, e and ig one can tell whether a pulse 
which arrives at the barrier is able to cross it. Moreover, by counting the 
number of maxima of activator within a certain time interval one can measure 
the frequency of the input and output signals (f, and f, respectively). When 
calculating fp and fọ a few initial pulses (usually 10 or 20) are neglected 
in order to eliminate a transient behavior at the beginning of evolution. 
To describe quantitatively the passive barrier as a device which transforms 
chemical signal frequency, the filtering ratio defined as f,/f, is introduced. It 
corresponds to the firing number, if notation of [65]-[69] is used. 

Calculations have been performed for different values of t, (or f,) and 
different values of the barrier’s width d. The diagrams in the space of param- 
eters (d,t,) showing ”phases” in which the barrier transforms the chemical 
signal in a given way are shown in figures presenting results obtained for 
particular models. 

The parameters d and t, (or f,) seem to be proper variables to describe 
how a passive barrier works as a transformer of the signal frequency. The 
interval of t, (or f,) is limited by the fact that once the excitable medium 
has been excited it needs some minimal amount of time (called the refractory 
period) in order to relax, before it may be excited again. Thus, there exists 
a minimal time pace tp min at which the pulses can be successfully initiated. 
The time t) min depends on the strength of excitation and it fixes the upper 
frequency of the signal. On the other hand, making t, very long reduces the 
problem to propagation of single pulses. The range of d is also finite as there 


always exists a barrier of width dmin, narrow enough to be transparent to 
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all pulses in a train and another one, impenetrable for any of them (of the 
width of daz). Thus, the nontrivial transforming properties of the passive 
barrier may be observed only in some finite range of £, and d, which depends 


on the selected model and values of its parameters. 


2.4.2 Results for the FH-N model. 


In calculations for the FH-N model the following parameters of a uniform 
grid have been used: n=251, n,=150, na=l160, 4=148, i2=162, 23=248. In 
order to produce a pulse the value of v has been decreased to vjn;=-0.8 on 
the left end of the interval. 

For the considered values of parameters the minimum time after which 
the second pulse may be re-excited is ty min,ru- N=2.4. However, for such a 
short time of consecutive excitations it is not possible to obtain a regular, 
stable train of pulses, because it happens that the medium at the initiation 
point is not well relaxed and an attempt to re-excite it by the assumed Vini 
fails. For the FitzHugh - Nagumo system the value of b min,FH-N Strongly 
depends on v;; (the strength of excitation) and also on dl, which describes 
the spatial size of excitation. Here in order to create a stable train of pulses 
tp€[2.90, 4.04] has been used . 

The maximum time within which the evolution is studied (tma.) is fixed 
in the calculations. The system is excited approximately p=tmar/tp times, 
so it is expected that such number of pulses is produced. The time shift 
between individual pulses within a regular train equals t, and f,=1/t, gives 
the frequency of incident pulses at the first indicator (4). 

To investigate the effect of the barrier tma,=500 and the time integration 
step dt,=5-10~% have been used to study a range of t, from 2.98 to 4.04 with 
increment 0.02 (f,e|0.248, 0.336]) and the barrier’s widths de[0.126, 0.172]. 
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For the selected tna, and tp, from 120 to 160 pulses have been observed in 
the system for each combination of tp and d. 

Figure 17 presents a sample signal (concentration of u) observed at the 
grid points 4, % and tz for dl=0.0175 (uniform grid) and f,=0.294. For all 
the results shown in Fig. 17 nj=150, while m is changed to obtain barri- 
ers of various width d. The upper curve (1) corresponds to incident pulses 
(reference signal at indicator i,). In case of m=151 (no passive barrier) the 
same signal is observed at indicator b with only a small time shift, which 
the signal needs to cover the distance between indicators 4 and iz (signal 
at ig - curve 2). The result is the same for na=155, which corresponds to a 
thin, fully transparent passive barrier (curve 3 presents the signal at i2). All 
incident pulses observed at indicator 7, get through the barrier and are also 
observed at indicator ię. However, for a properly chosen width the barrier 
becomes selective. For na=160 (the barrier’s width is then d=0.158) exactly 
every second of the incident pulses is transmitted through the barrier and 
observed at indicator b (curve 4) or indicator 43 (curve 5). 

Figure 18 summarizes the signal transforming properties of a passive bar- 
rier as a function of the barrier’s width d and the time shift between consec- 
utive incident pulses tp. The white, labeled regions correspond to the given 
filtering ratio (f,/f,). The area labeled as ”1” indicates the ratio equal to 
1, which means that every incident pulse is able to get through the passive 
barrier (the barrier is transparent to all pulses). When tp decreases - one 
arrives at the area (labeled as ”1/2” in Fig. 18) where only every second of 
the incident pulses is transmitted. With increasing d the filtering ratio de- 
creases, which means that the pulses are less and less frequently transmitted 
(see area ”1/3” in Fig. 18, where only one out of three of incident pulses 


gets through the barrier). Finally the barrier becomes too wide and no pulse 
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Figure 17. The time evolution of activator’s concentration at the grid point 
i (the upper curve 1), at the grid point ię (curves 2-4) and at the grid point 73 
far behind the barrier (curve 5). The evolutions at 72 correspond to no barrier 
(curve 2), a narrow transparent barrier (curve 3) and to a barrier for which 
frequency transformation occurs (curve 4). In the last case the evolution at 
iz is also shown (curve 5). Results for the FitzHugh - Nagumo model (Eqs. 
(2.1-2.4)), t,=3.40 (f,=0.294), d=0.070 (curve 3), d=0.158 (curves 4 and 5). 
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Figure 18. Filtering ratio (f,/f,) for the FitzHugh - Nagumo model as a 
function of the barrier’s width (d) and the time shift between consecutive 
pulses (tp). The white, labeled areas correspond to the situation when f, is 
the fraction of f, given in the picture. Gray color marks more complicated 
transformations of frequency. The dashed lines correspond to: (a) t,=3.10 
(fp=0.323), (b) t$=3.70 (f,=0.270), (c) d=0.160. The rectangular area in the 
bottom right hand side corner of the picture has been studied more carefully 


and the results are illustrated in Fig. 19. 
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can cross it (thus f,=0). This corresponds to the area labeled as ”0” in Fig. 
18.14 The dashed lines in Fig. 18 mark t,=3.10 (f,=0.323) - line (a), t,=3.70 
(fp=0.270) - line (b), and d=0.160 - line (c). The filtering ratios as functions 
of parameters from the dashed lines (a), (b) and (c) are discussed in Sections 
2.4.3 and 2.4.4. It is worth noticing that for the FitzHugh - Nagumo model 
the frequency transforming (corresponding to elimination of some of the in- 
coming pulses by the barrier) may appear for barriers which are narrower or 
wider than the penetration depth for a single pulse (de ru_N=0.163). In fact, 
in Fig. 18 there is a hidden line, corresponding to d=0.163, which divides the 
”1/2” area into two qualitatively different sub-areas. For d<0.163 the barrier 
is transparent for the first of incoming pulses, impenetrable for the second 
one etc. For d>0.163 the first of incoming pulses is stopped at the barrier, 
the second one gets through etc. In both cases the overall effect is the same 
(f,/fp=1/2). Thus, the described change occuring sharply for d=0.163 can 
be considered as phase shift of the outgoing signal. 

The part of the parameter space for which the filtering ratio is smaller 
than 1 has been studied more carefully with tmaz=500 and dt=2-10~%. The 
range of £, from 3.00 to 3.40 with increment of 0.01 and the range of d 
from 0.168 to 0.173 with increment of 0.00018 have been investigated (this 
is the rectangular area marked with dashed line in the bottom right hand 


side corner of Fig. 18). The results are presented in Figure 19, for which the 


Various modes of the frequency transformation occurring in the FH-N model are pre- 
sented on animated figures containing 6 incoming pulses each. The animations can be 
found on the CD-ROM (folder Figures). Fig_18a.htm presents the case of a fully trans- 
parent barrier, for which f,/f,=1. The width of the barrier is then increased. Fig_18b.htm 
shows fo/fp=2/3, Fig_18c.htm corresponds to f,/fp=1/2 and finally Fig_18d.htm presents 
an impenetrable barrier for which f,/f,=0. The frequency of the input signal f,=0.270 


(fp=3.70) is the same for all the figures mentioned here. 
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meaning of the colors and labels is the same as in Fig. 18. 

The gray regions in Figs. 18 and 19, which separate the areas labeled as 
"1", 1/2”, ”1/3”, correspond to more complex (periodic or non-periodic) 
transmission patterns described later in text. One of such patterns, corre- 
sponding to the filtering ratio f,/f,=2/5, is presented in Figure 20. Results 
presented in Fig. 20 have been obtained for n=251 (uniform grid), dl=0.01896, 
m=150 and m=160 (so the barrier’s width is d=0.1707). Indicators have 
been located at grid points i1=148, %=162 and i3=248. Curves 1, 2 and 3 
show the signal at indicators t, iz and iş respectively, for f,=0.318 (t,=3.14). 
Just after crossing the barrier the signal has the period of 5-t, with two max- 
ima of concentration per period (Fig. 20, curve 2). The signal far behind 
the barrier is shown as the curve 3. One can see that the small maxima 
have not developed into regular pulses and the signal is composed of peaks 
separated by 2-t, and 3-¢,. Such form of a signal has a transient character. 
The numerical experiments with pulses in the FH-N system have indicated 
that the velocity of a pulse increases with the distance separating it from the 
preceding ones [125]. Thus, a pulse which propagates 3-t, after the preceding 
one is faster than that propagating 2-t, after the preceding one. In the long 
time limit one obtains a periodic signal, the frequency of which is 2/5 of the 
original one. The stability of the 2/5 transforming mode described here has 


been confirmed by calculations carried out up to tmaz=10000. 


2.4.3 Numerical difficulties. 


The solutions which correspond to the filtering ratio 1/2 or 1/3 are quite 
stable numerically and they cover a large part of the parameter space. In 
order to learn more about what happens in the gray areas of Figs. 18 and 19 


a series of calculations has been performed with one of the parameters (the 
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Figure 19. Filtering ratio (f,/f,) for the FitzHugh - Nagumo model as a 
function of the barrier’s width (d) and the time shift between consecutive 
pulses (tp). The white, labeled areas correspond to the situation when f, is 
the fraction of fp given in the picture. Gray color marks more complicated 


transformations of frequency. 
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Figure 20. The time evolution of activator’s concentration in the FitzHugh 
- Nagumo model at the grid point 4 (curve 1), at the grid point 4 (curve 2), 
and at the grid point 73 (curve 3) for d=0.1707. Curves 1, 2 and 3 correspond 


to f„=0.318 and illustrate the frequency transforming ratio equal to 2/5. 
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width of the passive barrier d or the time shift between consecutive incident 
pulses t,) fixed and the other changing within a certain range. By monitoring 
the input and output signals (similar to those presented in Fig. 17) one is 
not only able to calculate the filtering ratio for the given combination of d 
and t,, but also describe the sequence of transmitted/stopped pulses in the 
input signal, which leads to a particular value of f,/f,. The calculations have 
been done for parameters on the lines of constant d or t, shown in Fig. 18. 
The computations have to be performed for a period of time long enough 
to observe many full cycles of signal transformation. However, it has been 
discovered that the results of calculations are quite sensitive with respect to 
the size of the grid and the time integration step (dt) used in computations. 
Figure 21 presents a set of curves showing the filtering ratio for t,=3.10 
(fp=0.323) as a function of d. The considered values of parameters line on 
the horizontal dashed line (a) in Fig. 18. All results have been obtained 
for dt=1-107% and the calculations have been performed up to tmad=500. 
For curve (a), marked with empty triangles, a uniform grid with n=400 has 
been used. Curve (b), marked with empty diamonds, indicates the results 
for a uniform grid with n=800. The passive gap has been located between 
ny =180 and m=189 for n=400, n,=360 and m=377 for n=800. In case of 
curves (c) - filled circles and (d) - empty circles, adaptive grids with n=800 
(c) or n=870 (d) respectively have been applied. The value of dł in these 
calculations depends on system’s length l and the number of grid points 
used. The signals at indicators 4 and tz, located at approximately constant 
distance of l;=0.6 (in the dimensionless units of distance) before and behind 
the passive barrier have been analyzed (cf. the bottom part of Fig. 15). 
Results presented in Fig. 21 show a rich structure of different filtering 


ratios, which is much more complex than the one shown in Fig. OC for 
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Figure 21. Results for the FitzHugh - Nagumo model. Filtering ratio fo/fp for 
fixed t,=3.10 (/,=0.323) and d changing within the presented range. All results 
have been obtained for dt=1-107*, tmar=500, but for each of the curves (a) - (d) 
different grids have been used. This figure illustrates the dependence of the results 
on computational parameters. Curve (a) - marked with empty triangles - uniform 
grid with n=400. Curve (b) - marked with empty diamonds - uniform grid with 
n=800. Curve (c) - marked with filled circles - adaptive grid with n=800. Curve 


(d) - marked with empty circles - adaptive grid with n=870. 
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the ” naive” model of an excitable system. Of course, one can observe some 
similarities: the region of d for which f,/f,=1/2 is the dominant one, and 
the second most important is f,/f,=1/3. But there are also filtering ratios 
between 1/3 and 1/2, which are absent in the ”naive” model (see Section 
1.2.2). The calculations show that although f,/f, as a function of d looks 
very similar for different values of parameters of integration, it shifts towards 
greater values of d when more accurate integration techniques are applied. A 
significant shift between curves (a) and (b) in Fig. 21 indicates that for the 
uniform grids the results strongly depend on the grid size. It is worth noticing 
that the behavior of the filtering ratio as a function of d does not change, but 
the function as a whole is just shifted towards wider barriers. It is expected 
that for yet finer grid one should obtain results which are numerically stable, 
but a finer grid means that more grid points should be used to describe a 
system of the same size. However, in the implicit method of solving parabolic 
reaction-diffusion equations the solution at each time step is obtained via 
iterations. In case of the FH-N model the roundup errors lead to instabilities 
when few thousands of grid points are used. The adaptive grids (Fig. 16) 
allow one to obtain more accurate numerical solutions. For such grids the 
fine resolution within the most important area of the investigated system may 
be achieved without using large number of grid points, which provides both 
accuracy and stability of computations. As already mentioned, curve (c) in 
Fig. 21 has been obtained for the adaptive grid with n=800. The resolution 
in the neighborhood of the passive gap for this grid is dl,=dl/8~0.06 (for 
1=25.72), which is four times smaller than dl for the finest uniform grid. For 
curve (d) this resolution is twice higher, but f,/f,(d) obtained for both grids 
(curves (c) and (d) in Fig. 21) are almost identical. Therefore, the adaptive 


grid with n=800 seems to be sufficient for the calculations and it has been 
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used to obtain results given below. 

Several values of the time step of integration (df=5-10~%, 1-1073 and 
1-1074) have been used to verify the consistency of the results. Although the 
implicit algorithms should be in general stable for all values of time and space 
integration steps, it has been observed that numerical instabilities may ap- 
pear for dt=5-10~? and large n. On the other hand, the results for dt=1-107° 
and dt=1-107* have been regular and consistent. Consequently, dt=1-10~3 
has been used to study the complex frequency transforming patterns, de- 


scribed in the next section. 


2.4.4 Complex patterns of output signals. 


The results presented in this section have been obtained for the adaptive grid 
with n=800. The passive area has been located between n,=400 and ng=431 
(30 grid points inside the passive barrier). The incoming and outgoing pulses 
have been monitored at grid points 7;=390 and %=441 respectively. The 
computations have been carried out with the time step dt=1:107*% up to 
tmax=000, if not explicitly stated otherwise. Circa 150 pulses arrive at the 
barrier within the selected time tmaz- 

A special notation has been introduced to describe the output signal 
[124]. The incident pulses observed at 4 introduce a natural time scale in 
the system. One writes ”1” if a pulse gets successfully through the passive 
barrier and is observed at 72 or ”0” otherwise. Thus the output signal may be 
coded as a sequence of ”0” and” 1”. In such notation a common case in which 
every second pulse passes (filtering ratio f,/f,=1/2) is coded as (01) and the 
mode ”1/3”, is described as (001). It is understood that the given sequence 
repeats periodically. A pattern coded as (abc...)P(def...)? (where p and q are 


positive integer numbers) means, that behind the barrier first the sequence 
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(abc...) is observed p times, and next the sequence (def...) appears q times. 
In this notation non-periodic modes correspond to an infinite sequence of 
(abcde...). If pattern’s sequence is finite, then it is very easy to calculate the 
corresponding filtering ratio f,/f, because it equals to the sum of the symbols 
in the sequence divided by the number of symbols. For non-periodic modes 
one can estimate the filtering ratio using a finite part of the sequence and 
of course the more pulses is taken into account, the better approximation is 
achieved. 

Figure 22 presents the filtering ratio for t,=3.10 (f,=0.323) as a function 
of the width of the passive barrier changing from 0.171 to 0.176 with incre- 
ment of 0.0001 (this has been achieved by changing the total length of the 
system | from 27.49 to 28.19 with increment of 0.016). This value of t, corre- 
sponds to line (a) in Fig. 18. The results coming from our computations are 
marked with filled circles. In Fig. 22 one may notice several plateaus, which 
correspond to simple patterns of frequency transforming and are labeled with 
appropriate filtering ratios. Looking from left to right there are modes: ” 1/2” 
with the sequence (01), ”2/5” with (01)(001), ”1/3” - (001), ”1/4” - (0001) 
and finally ”0” with the sequence (0). It is remarkable that a ”simple” fre- 
quency transformation occurs within a wider range of barrier’s width than 
complex ones. Other interesting transmission patterns, located between those 
plateaus, have also been observed, but the ranges of values of d within which 
those patterns appear are so narrow that they have been detected just for 
a single width of the barrier. The unique patterns have been labeled with 
letters a - k in Fig. 22. Probably the most complicated transmission pat- 
tern is associated with point a, located just below the plateau ”1/2”. In 
this case the sequence of pulses behind the barrier is: (01)'*(001)(01)" (001) 
and the corresponding filtering ratio f,/f,=29/60. For point b in Fig. 22 
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Figure 22. Results for the FitzHugh - Nagumo model. Filtering ratio fo/fp for 
fixed tp=3.10 (f,=0.323) and d changing within the presented range. The filled 
circles correspond to computational points. The curve forms several plateaus, 
which are labeled with corresponding values of the filtering ratio (1/2, 2/5, 1/8, 
1/4 and 0). At points labeled a - k the following values of f,/f, (and trans- 
mission patterns) have been observed: a - 29/60, (01)'*(001)(01)'8 (001); b - 
11/24, (01)4(001)(01)° (001); c - 4/9, (01)? (001); d - 3/7, (01)? (001); e - 9/22, 
(01001) (0101001); f - 3/8, (01)(001)?; g - 7/20, (01)(001)$; h - 2/7, (001)(0001); 
j-1/5, (00001) and k - 1/7, (0000001). 
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the corresponding pattern is (01)*(001)(01) (001) resulting in f,/f,=11/24, 
point c - (01)? (001), fo/fo=4/9; point d - (01)? (001), f,/f,=3/7; point e - 
(01001)? (0101001), f,/f,=9/22; point f - (01)(001), fo/fp=3/8; point g - 
(01)(001)°, fo/fp=7/20; point h - (001)(0001), f/f,=2/7; point j - (00001), 
fo/fp=1/5 and point k - (0000001), f,/f=1/T. 

Figure 18 suggests that transmission patterns corresponding to filtering 
ratios greater than 1/2 are also present. Such patterns are absent in the 
” naive” model of excitable systems (Figs. OB - 0C). To see them clearly 
the calculations for t,=3.70 (f,=0.270) have been performed for the barrier's 
width d changing from 0.130 to 0.171 with increment of 0.001 (l changing 
from 20.93 to 27.49 with increment of 0.16). This value of t, corresponds to 
line (b) in Fig. 18. The results are plotted in Figure 23. Indeed, one may 
see here several plateaus, labeled with corresponding filtering ratios, which 
are greater than 1/2. From left to right there are: plateau ”1” with the 
sequence (1); plateau ”6/7” - (0111111); plateau ”4/5” - (01111); plateau 
»3/4 - (0111); plateau 2/3” - (011); plateau ”1/2” - (01) and plateau 
”0” with the trivial transmission pattern (0). Even more complex output 
signals have been observed. The points labeled with letters correspond to 
the following patterns and filtering ratios: point a - (0)(1)*, f,/f,=35/36; 
point b - (0)(1)™, fo/f»=14/15; point c - (0)(1)", fo//,=10/11; point d - 
(0)(1), fo/fo=8/9 and point e - (011111), f,/f(,=5/6. The output signals 
observed here are dual to those with filtering ratios smaller than 1/2. 

Figure 24 presents the filtering ratio for d=0.160 (1=25.72) plotted versus 
the frequency of incident pulses f,e|0.25, 0.28] (the time shift between consec- 
utive pulses changing from 3.57 to 4.00 with increment of 0.02; approximate 
values of t, are given on the top axis of Fig. 24). These values of d are placed 


on the line (c) in Fig. 18. The numerical labels in Fig. 24 give the filtering 
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Figure 23. Results for the FitzHugh - Nagumo model. Filtering ratio fo/fp 
for fixed t,=3.70 (f,=0.270) and d changing within the presented range. The 
filled circles correspond to computational points. The plateaus are labeled with 
corresponding values of the filtering ratio (1, 6/7, 4/5, 3/4, 2/3, 1/2 and 0). At 
points labeled a - e the following values of fo/fp (and transmission patterns) have 
been observed: a - 35/36, (0)(1)%; b - 14/15, (0)(1)“; c - 10/11, (0)(1)'°; d - 
8/9, (0)(1)$ and e - 5/6, (011111). 


CHAPTER 2. RESULTS. 85 


Approx. time shift between consecutive pulses (tp) 


4.00 3.90 3.80 3.70 3.60 
| | | | | 


Ratio of frequencies (f/ fp) 


d=0.160 D 


0.25 0.26 0.27 0.28 
Frequency of incident pulses (fp) 


Figure 24. Results for the FitzHugh - Nagumo model. Filtering ratio fo/fp 
for fixed d=0.160 and fp changing from 0.25 to 0.28 (see the bottom axis), which 
corresponds to tp decreasing from 4.00 to 3.57 (see the approximate top axis). The 
filled circles correspond to computational points. The plateaus as well as single 
points are labeled with corresponding values of fo/fp (1, 9/10, 4/5, 3/4, 2/3 and 
1/2). The letters A, B, C and D mark the intervals of the diagram which have been 
studied with a higher resolution. The results for the interval AB are presented in 


Fig. 25 and for the interval CD in Fig. 26. 
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ratios observed for corresponding plateaus or points. In Fig. 24 from left to 
right there are: plateau ”1” with the sequence (1); point ”9/10” with the 
sequence (0)(1)?; point ”4/5” - (01111); point ”3/4” - (0111); plateau ” 2/3” 
- (011) and plateau ”1/2” with (01). In order to check whether the filtering 
ratio is a monotonic function of fp, a number of calculations for f„e|[0.2604, 
0.2618] and for f,e|0.2646, 0.2674] have been performed. The ends of these 
intervals correspond to fp for points A, B, C, D in Fig. 24. The system’s 
behavior has been probed using 6t,=0.002 (corresponding to 6f,=0.00014). 
In the interval AB the behavior is trivial, which means that studied systems 
splits into two classes corresponding to (0111) and (011) modes (cf. Figure 
25). 

Similar calculations have been performed for systems located between 
points C and D in Fig. 24 (this interval corresponds to f,e|[0.2646, 0.2674]). 
In this case two patterns (a and 6 in Figure 26) have been observed between 
the plateaus corresponding to filtering ratios 2/3 ((011) mode) and 1/2 ((01) 
mode). One of these patterns (5) describes quite simple transformation of the 
original signal to (01)(011) mode, which is just 1:1 mixture of the neighboring 
modes. The other point (a) corresponds to more interesting behavior. Using 
the same dt and dl, as for the other points shown in Fig. 26 (filled circles) 
the filtering ratio f,/f,=0.51917 has been observed. The decrease in time 
step dt increases the filtering ratio as follows: dt=1-10~ gives f,/f,=0.55602 
(empty circle), dt=2-10~° gives f,/f,=0.57594 (empty triangle), dt=1-107* 
gives f,/fp=0.57895 (cross). Moreover it has been observed that time intervals 
between transmitted pulses are not always a multiplicity of t,. In order to 
explain it the time evolution u(t) at the points 4 and p has been examined. 
Figure 27 shows the relevant part of it. One can notice that in some cases 


the transmitted signal does not develop into a pulse, u(t) slightly decreases 
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Figure 25. Results of a more detailed examination of the interval AB from Fig. 
24. Filtering ratio fo/fp for fixed d=0.160 and fp changing from 0.2604 to 0.2618. 
Points A and B in this figure correspond exactly to points A and B in Fig. 24. 


No filtering ratios different from 3/4 or 2/3 have been observed. 
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Figure 26. Results of a more detailed examination of the interval CD from Fig. 
24. Filtering ratio f,/f, for fixed d=0.160 and fp changing from 0.2646 to 0.2674. 
The filled circles mark the computational points. Points C and D in this figure 
correspond exactly to points C and D in Fig. 24. A "strange" point a has been 
found between plateaus 2/3 and 1/2 (labeled with values of the filtering ratio), 
for which the value of f/f, depends on the time step of integration. Different 
symbols mark the values of the filtering ratio obtained for: dt=1-10~° (filled circle), 
dt=1-107* (empty circle), dt=2-10~° (empty triangle), dt=1-10~° (cross). In the 
studied region the point @ has been observed, which is regular and corresponds to 


the transmission pattern (01)(011) (thus f,/f,=3/5). 
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Figure 27. Results for the FitzHugh - Nagumo model. A part of the time 
evolution of the value of u at the grid point i (input signal - the upper curve 1) 
and at the grid point i (output signal - curve 2) observed for d=0.160, fp=0.2660 
with dt=1-10~° (point a - the cross in Fig. 26). The arrow marks the position of 


the ”strange excitation” in the output signal. 
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and then it starts to increase again, forming a pulse shifted with respect 
to the forcing signal. Such strange behavior for f,=0.2660 (t,=3.76) does 
not disappear when decreasing dt and dl, The moment it appears first 
and the time intervals between successive strange excitations depend on dt, 
but no regularities have been noticed. The presence of irregular ”strange 
excitations” makes the filtering ratio at the point a dependent on dt. The 
value of f,/f, at this point still remains smaller than 0.6, so one cannot exclude 
a non-monotonic dependence of the filtering ratio on f,/f,. Moreover, such 
strange excitations have been observed for the whole range of time integration 
steps used, so one is unable to blame numerical instabilities for their presence. 
The time evolution of the output signal between the strange excitations is a 


mixture of (01) and (011) modes with different proportions. 


2.4.5 Results for the R-Z model. 


In the calculations for the Rovinsky - Zhabotinsky model a uniform grid 
with n=320, ny=150, no=155 has been used as well as the dimensionless 
time integration step dr=1:10%. The pulses have been initiated at the left 
end of the interval by increasing the value of z to 0.1. 

Figure 28 presents a typical signal (value of x) observed at the first and 
second indicators (grid points 4 =148 and b=157) for dl=0.814 and 7,=30.0 
(255sec). For these values of parameters the excitation at the boundary gives 
a regular and stable train of pulses. In calculations, the results of which are 
shown in Fig. 28, m is changed to obtain barriers of different widths. The 
upper curve (1) corresponds to incident pulses (reference signal at indicator 
1). For ng=151 (no passive barrier - curve 2) the same signal (shifted in time) 
is observed at indicator 2. The same behavior is observed for a thin, fully 


transparent passive barrier (n2=153 - curve 3). All incident pulses observed 
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Figure 28. The time evolution of activator’s concentration at the grid point 
i (the upper curve 1) and at the grid point % (curves 2-4). The evolutions at 
i2 correspond to no barrier (curve 2), a narrow transparent barrier (curve 3) 
and a barrier for which frequency transformation occurs (curve 4). Results for 
the Rovinsky - Zhabotinsky model (Eqs. (2.5-2.8)), Tp=30 (fp=0.033)), d=1.628 
(0.0150- Dzem) (curve 3), d=3.256 (0.0300- Dzem) (curve 4). Tp=30 corre- 


sponds to 255sec. 
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at indicator 1 get through the barrier and are also observed at indicator 2. 
For no=155 every third of the incident pulses is transmitted through the 
barrier as observed at indicator 2 (curve 4). In this case the width of the 
passive barrier is 3.256 (0.0300- DZ ccm). 

It has been found that for the > model the minimal time necessary to 
initiate a new pulse after the first pulse has been produced in the system 
IS Tpmin,R-z%5.8 (49sec) but one has to use much higher values of 7, to 
obtain a stable, regular train of pulses. The range of r, from 30 (255sec) 
to 60 (510sec) with increment of 2 (17sec) and the range of d changing 
from 3.188 (0.02939- Dzem) to 3.295 (0.03037- Dzem) have been stud- 
ied. Tmax 2000 (17000sec) has been used, so during the evolution from 33 
to 66 pulses may appear in the system for each combination of T, and d. 

The results are shown in Figure 29. For all the examined values of d and 
Tp every excitation creates a pulse. The filtering ratio f,/f, presented in Fig. 
29 reveals areas of parameters’ values for which every second or every third 
of the incident pulses is transmitted (these regions are labeled with ”1/2” 
and ”1/3” respectively). Label ”0” in Fig. 29 indicates the area in which 
no pulse can cross the passive barrier. The gray regions between the white, 
labeled areas in Fig. 29 correspond to more complex transmission patterns. 
Unlike for the FH-N model, the transformation of frequency occurs only for 
barriers which are narrower than the penetration depth for a single pulse 
given in Section 2.2 d, p_7%3.295 (0.03037-, / #5 xem). 

Figure 30 presents the filtering ratio fo/ fp a versus the frequency of 
incident pulses f, for a selected barrier's width d=3.2375 (0.0298- Dzem). 
This value of d corresponds to the dashed vertical line in Fig. 29. Here the 
frequencies f, and f, are dimensionless frequencies, calculated as the inverse 


of dimensionless time r. For f,€|0.0167, 0.0278] every second of incident 
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Figure 29. Filtering ratio (f,/f,) for the Rovinsky - Zhabotinsky model as 
a function of the barrier’s width (dimensionless - d) and the time shift between 
consecutive pulses (dimensionless - Tp). The white, labeled areas correspond to the 
situation when fo is the fraction of fp given in the picture. Gray color marks more 


complicated transformations of frequency. The dashed line indicates d=3.2375. 
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Figure 30. Filtering ratio in the Rovinsky - Zhabotinsky model for a selected 
barrier's width d=3.2375 (0.0298- Dzem). This value of d corresponds to the 
dashed vertical line in Fig. 29. f,/fp is presented as a function of the dimensionless 
frequency of incident pulses (f, - bottom axis) or the (approximate) physical time 
shift between consecutive incident pulses (tp [sec] - top axis). Labels give the 


filtering ratio (fo/ fp). 
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pulses gets through the barrier. For f,>0.0294 only one out of three incident 
pulses gets through the barrier (thus the filtering ratio is f,/f,=1/3). 


2.4.6 Results for the Oregonator model. 


In order to check whether the frequency transforming is a typical feature of a 
passive barrier in an excitable medium calculations for the Oregonator model 
[126, 127, 128] have been performed. The form of the model presented in [127] 
with diffusion of the activator added has been used. In the Oregonator model 
the active field is described by the following reaction - diffusion equations 


(2.17-2.18): 


A 2 |od bf aja eea 2 
7: zul u)— f-v TE + DV*u (2.17) 
m =u-v (2.18) 


where u corresponds to the scaled concentration of activator (HBrO>) and v 
- to the scaled concentration of catalyst (Ce**) [127]. In the passive regions, 
without catalyst, the concentrations of u and v evolve according to Eqs. 
(2.19-2.20): 

"u" + DV?u (2.19) 


= const. (2.20) 


In Eqs. (2.17-2.20) t denotes the scaled time and e (a time scale parameter) 
is small [127]. Another small parameter q is connected to the rate constants 
of the reactions involved in the Oregonator model. The stoichiometric pa- 
rameter f is proportional to the average number of bromide ions released per 
metal ion reduced by organic matter. The details of scaling are described in 
[127]. 

In the calculations the following values of parameters have been used: 


e=0.05 [128], f=3 [127], q=0.0002 [127] and D=1.0. For these values of 
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parameters the stationary state of the system within the active areas corre- 
sponds to 


Usa = Usa = 3.9988 - 1074 (2.21) 


and it is excitable. In the passive areas the stationary solution is given by 
Us = t= 0 (2.22) 


In the calculations for this model a uniform grid with dl=0.3733, (n=400), 
m=200, dt=1-10~* and t,=10.4. The pulses have been initiated on the left 
end of the interval by decreasing v to vjn;=0.0. The evolution of this system 
has been studied up to tmar=1000 (so that over 95 pulses has been involved 
in each experiment). It has been observed that a passive barrier in the 
Oregonator model also reveals frequency transforming properties. Figure 
31 presents a sample signal (the value of u) observed at the first and second 
indicators (grid points 4=198 and 4=213 respectively). As mentioned above, 
m=200 while n is changed to obtain barriers of different width. The upper 
curve (1) corresponds to incident pulses (reference signal at indicator 1). 
For n2=201 (no passive barrier - curve 2) the same signal (shifted in time) 
is observed at indicator 2. The same behavior is observed for a narrow, 
fully transparent passive barrier (ng=206 - curve 3). All incident pulses 
observed at indicator 7, get through the barrier and are also observed at 
indicator 4. For a wider barrier (nz=211) every second of the arriving pulses 
is transmitted through the barrier and may be observed at indicator 2 (curve 
4). For curve 4 the width of the passive barrier is 3.733. In this case, as for 
R-Z (BZ) system, the frequency transforming occurs for barriers narrower 
than the penetration depth for a single pulse, because the first transmitted 
pulse makes the barrier impenetrable for the subsequent ones (note that in 


Fig. 31 the first of incident pulses gets through the barrier). 
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Figure 31. The time evolution of activator’s concentration at the grid point 4 
(the upper curve 1) and at the grid point ię (curves 2-4). The evolutions at b 
correspond to no barrier (curve 2), a narrow transparent barrier (curve 3) and a 
barrier for which frequency transformation occurs (curve 4). Results for the Oreg- 
onator model (Eqs. (2.17-2.20)), tp=10.4 (fp=0.096), d=1.867 (curve 3), d=3.733 


(curve 4). 
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2.4.7 Implications of frequency transforming for se- 


lected signal processing reactors. 


Having in mind the results of previous sections one may expect that the 
properties of signal processing reactors which include a passive barrier should 
depend on input signal frequency. In this chapter the effect of frequency of 
input signal on two of such devices is discussed. 

A chemical diode is an asymmetric junction of two active areas, which 
transmits pulses of excitation coming from one direction and stops those ar- 
riving from the other one. The structure of the diode, as proposed by I. 
Motoike and K. Yoshikawa in [86] is shown schematically in Figure 32, where 
the darker area stands for the passive medium and the brighter areas cor- 
respond to the active one. The left edge of the passive region is a straight 
line, while the lines bounding of the active region on the right are perpen- 
dicular and form a tip. The shortest distance between the active areas is 
d, >0. No flux boundary conditions at all borders of the structure presented 
in Fig. 32 and free flow of mobile reagents between active and passive areas 
are assumed. Plain pulses of excitation may be initiated at the left hand side 
border and travel to the right or they can be excited on the right hand side 
border and travel to the left. The dashed lines on both sides of the junction, 
denoted by 7, and ip, mark the indicators. The propagation of pulses through 
the junction is studied by measuring the mean value of the concentration of 
the activator (u or x) along these lines. Indicators i, and ip are located 
symmetrically with respect to the longitudinal axis of the junction and their 
length is about half of the width of the whole device. A travelling pulse of 
excitation may get from one active area to the other one through the passive 
gap, provided that the gap is not too wide. 

As already mentioned (Section 2.2), I. Motoike and K. Yoshikawa found 
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Figure 32. The scheme of a chemical diode. Brighter areas correspond to active 
regions, darker ones show the passive gap. The gap starts to the right from line 
nı of the grid (excluding nı). Two dashed lines on both sides of the junction show 
the position of ” indicators”, at which the mean values of u (for the FH-N model) 
or z (for the R-Z model) as functions of time are recorded. For such geometry of 
the diode a single pulse of excitation may travel through the junction from left to 


right, but not in the opposite direction. 


CHAPTER 2. RESULTS. 100 


[86] that the penetration depth strongly depends on the geometry of the 
boundary between active and passive media. They used the FH-N model 
(Eqs. (2.1-2.4)) and studied the propagation of pulses in a junction shown 
in Fig. 32. According to their results [86], if the penetration depth for 
a pulse propagating to the right is denoted by de, then the penetration 
depth for a pulse travelling in the opposite direction is only 0.63:d,. Re- 
sults obtained for the passive barrier within the FH-N model (Section 2.2) 
give de py_n¥0.163. Now building a diode is very easy, because it is sufficient 
to separate the active areas by a passive field characterized by appropriate 
dn, such that 0.63:d, pFH-N<dn,FH-N <de,rpu—n. For the purpose of tis work 
dn. FH-N=0.1152 has been selected. A similar geometry of the reactor, which 
works as a signal diode, may be used in the case of Rovinsky - Zhabotinsky 
model. Here dn R-z=3.0 (0.02765- Dzem) has been used. This is slightly 
smaller than the penetration depth d. R_z=3.295 (0.03037- pam). 

The work of a chemical diode is illustrated in Figure 33 (for the FH-N 
model) and Figure 34 (for the R-Z model). In both figures two upper curves 
(no. 1 and 2) correspond to a single plain pulse, excited on the left hand 
side boundary of the system and travelling to the right. The mean values of 
activator (u or z) along 4, (curve 1) and 4 (curve 2) are plotted as functions 
of time given at the top axis. The pulse which arrives first at 4, gets through 
the passive gap and a moment later it is observed at % (curve 2). The two 
lower curves (no. 3 and 4) illustrate the propagation of a pulse initiated on 
the right hand side boundary and travelling to the left. The mean values 
of u (or x) along b (curve 3) and 4 (curve 4) are plotted as functions of 
time shown on the bottom axis. In this case the pulse is first observed at ig 
(curve 3), but then it ”dies out” within the passive area (curve 4 is just a 


straight line). The arrows mark the direction of propagation of the pulses, 
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Figure 33. The time evolution of the value of u at the indicators 7, (curves 1 
and 4) and ig (curves 2 and 3). The two upper curves (no. 1 and 2) correspond 
to a pulse propagating from left to right (orientation defined by Fig. 32) and the 
corresponding time scale is at the top axis. In this case the diode transmits the 
pulse. The two lower curves (no. 3 and 4) correspond to a pulse propagating 
from right to left (cf. Fig. 32). In this case the pulse is stopped. The arrows 
mark the direction of propagation. Results for the FitzHugh - Nagumo model, 
dn, pH-N=0.1152. 


CHAPTER 2. RESULTS. 102 


Time (dimensionless) >>> 
0.00 20.00 40.00 


concentration of x 
V 
V 
V 
Po 


Signal - 


<<< 


| | 
40.00 20.00 0.00 
<<< Time (dimensionless) 


Figure 34. The time evolution of activator’s concentration at the indicators 4 
(curves 1 and 4) and % (curves 2 and 3). The two upper curves (no. 1 and 2) 
correspond to a pulse propagating from left to right (orientation defined by Fig. 
32) and the corresponding time scale is at the top axis. In this case the diode 
transmits the pulse. The two lower curves (no. 3 and 4) correspond to a pulse 
propagating from right to left (cf. Fig. 32). In this case the pulse is stopped. The 
arrows mark the direction of propagation. Results for the Rovinsky - Zhabotinsky 


model, dn. R-z=3.0 (0.02765: poem). 
0 
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where the orientation is defined by Fig. 32. In Fig. 34 the mean values of x 
at the indicators are plotted versus dimensionless time r (T=1 corresponds 
to 8.5sec)."$ 

The results presented in Figs. 33 - 34 have been calculated using the im- 
plicit method already mentioned. For the FH-N model a rectangular area of 
4.80x0.9792 covered with a square grid of 251x52 points has been considered 
(including borders, thus dl=0.0192). Indexing the vertical lines of the grid 
from 0 (left hand side border) to 250 (right hand side one) it may be speci- 
fied that indicators are placed along 44 =80 and p=126. Both of them include 
25 nodes. The passive area starts to the right from line mm =100 (this line 
of the grid is still characterized by the active dynamics). The shortest dis- 
tance between active areas is dp pu_N=0.1152, which is covered by six nodes. 
Time integration step dt=5-10~3 has been used. For the R-Z model the area 
of 240x38.25 has been covered with a square grid of 321x52 points, hence 
dl=0.75. The other parameters are respectively: 1=130, b=175, m=150 
and dn. R-z=3.0 (0.02765- Dzem - four nodes included). In this case the 
time integration step dt=1-107°. 

For the same values of parameters, calculations for a regular train of 
pulses arriving at a diode have been performed. The results are presented in 
Figures 35 - 37, which have exactly the same structure and symbols as Figs. 
33 - 34. 

For the FH-N model the input signal with frequency f,=0.278 (corre- 


sponding to the time shift between consecutive pulses t,=3.6) has been cho- 


'5The animated illustrations to Figures 33 - 34, showing the propagation of single pulses 
through the chemical diode, can be found on the CD-ROM (folder Figures). Fig_33a.htm - 
FH-N model, pulse propagating to the right; Fig_33b.htm - FH-N model, pulse propagating 
to the left (stopped by the diode); Fig_34a.htm - R-Z model, pulse propagating to the right; 
Fig_34b.htm - R-Z model, pulse propagating to the left (stopped by the diode). 
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sen. The time evolution of the values of u at 4 and is are plotted in Fig. 
35. As expected, the diode is transparent to all pulses travelling to the right 
(in the forward direction) - curves 1 and 2. However, for the chosen values 
of parameters the junction becomes transparent to some of the pulses trav- 
elling in the reverse direction (to the left - curves 3 and 4) . Therefore, the 
diode works as a transformer of chemical signal frequency.'* It is easy to 
explain this observation. In this case the barrier (dy ru-N=0.1152) is nar- 
rower than the penetration depth in the forward direction (d. rH-N=0.163), 
therefore each pulse coming from this direction is transmitted. On the other 
hand it has been selected as being only slightly wider than the penetration 
depth in the reverse direction: 0.63-d. pu_N=0.1027. Consequently, the gap 
demonstrates the same type of frequency transforming as the passive barrier 
discussed previously (see Section 2.4.2). 

For the R-Z model the input signal with dimensionless frequency f,=0.023 
has been used (this corresponds to the time shift between consecutive pulses 
Tp=44, ie. 374sec). The results are plotted in Figure 36. Comparing curves 
3 and 4 (in Fig. 36) one can see that in this case the diode stops all the 
pulses propagating in the reverse direction. As shown in Fig. 34, dn nR-z=3.0 
is larger than the penetration depth in this direction. The results of Section 
2.4.5 for the R-Z model (Fig. 29) indicate that such system becomes less 
” transparent” at higher signal frequencies. Therefore, it is natural that all 
pulses propagating in the reverse direction are stopped. However, the diode 
fails to transmit some pulses propagating in the forward direction (here: the 


2nd and 4th ones, see curves 1 and 2). This result can be easily explained 


‘©The animated illustrations to Figure 35, showing the propagation of this short train 
of pulses through the chemical diode, can be found on the CD-ROM (folder Figures). 
Fig_35a.htm - FH-N model, pulses propagating to the right; Fig_85b.htm - FH-N model, 
pulses propagating to the left (note that not all of them are stopped by the diode). 
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Figure 35. The time evolution of the value of u at the indicators i (curves 1 
and 4) and ig (curves 2 and 3). The two upper curves (no. 1 and 2) correspond 
to a train of pulses propagating from left to right (orientation defined by Fig. 
32) and the corresponding time scale is at the top. The two lower curves (no. 
3 and 4) correspond to a train of pulses propagating from right to left (cf. Fig. 
32). The arrows mark the direction of propagation. Results for the FitzHugh - 
Nagumo model, dn FH-N=0.1152, t,=3.60 (f,=0.278). Although the diode works 
still properly in the forward direction (curves 1 and 2), some pulses may cross it 


also in the reverse direction (curves,3 and 4). 
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Figure 36. The time evolution of activator’s concentration at the indicators 4 
(curves 1 and 4) and % (curves 2 and 3). The two upper curves (no. 1 and 2) 
correspond to a train of pulses propagating from left to right (orientation defined 
by Fig. 32) and the corresponding time scale is at the top. The two lower curves 
(no. 3 and 4) correspond to a train of pulses propagating from right to left (cf. 
Fig. 32). The arrows mark the direction of propagation. Results for the Rovinsky 
- Zhabotinsky model, dh,R-z=3.0 (0.02765- DCM), Tp=44 (f,=0.023). Tp=44 
corresponds to 374sec. In this case the diode stops every second pulse travelling 
in the forward direction (curves 1 and,2),. while it works correctly in the reverse 


direction (curves 3 and 4). 
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by what is already known from Section 2.4.5. Although dz R-z <de r-z (as 
3.0<3.295), the gap is still in the region for which the frequency transforming 
in the R-Z model occurs."” 

When the frequency of the input signal is increased - only the first pulse 
of the train is able to get through the diode in the forward direction and 
all subsequent pulses are stopped. This situation is shown in Figure 37, 
where the frequency of the incoming signal f,=0.05 has been used (7,=20, 
corresponding to 170sec)."* 

The response of the cross junction (Section 2.3, Fig. 3) to a train of pulses 
travelling inside a channel has also been studied. It has been found that if the 
time shift between consecutive pulses in a train is small, the characteristics 
of the junction changes (with respect to the results given in Sections 2.3.1 - 
2.3.2, concerning single pulses). A part of the cross junction has been studied, 
for the FH-N model with the inner width of the channel equal to 1.02 and 
the width of the passive stripes of 0.16 (exactly as in Section 2.3.1), using 
the implicit method with dl=0.02 and dt=5:10 %. The pulses inside one of 
the signal channels have been initiated with the time shift equal to t,=3.6. 
The result of calculations is presented in Figure 38, where brighter colors 


correspond to higher values of u. The first pulse behaves exactly as expected: 


it remains within the channel (as the width of the passive stripes is just above 


"The animated illustrations to Figure 36, showing the propagation of this short train 
of pulses through the chemical diode, can be found on the CD-ROM (folder Figures). 
Fig_36a.htm - R-Z model, pulses propagating to the right (note that the second and fourth 
of them are stopped by the diode); Fig_36b.htm - FH-N model, pulses propagating to the 


left. 
18The animated illustrations to Figure 37, showing the propagation of this short train 


of pulses through the chemical diode, can be found on the CD-ROM (folder Figures). 
Fig_37a.htm - R-Z model, pulses propagating to the right (note that all but the first one 
are stopped by the diode); Fig_37b.htm - FH-N model, pulses propagating to the left. 
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Figure 37. The time evolution of activator’s concentration at the indicators 4 
(curves 1 and 4) and % (curves 2 and 3). The two upper curves (no. 1 and 2) 
correspond to a train of pulses propagating from left to right (orientation defined 
by Fig. 32) and the corresponding time scale is at the top. The two lower curves 
(no. 3 and 4) correspond to a train of pulses propagating from right to left (cf. 
Fig. 32). The arrows mark the direction of propagation. Results for the Rovinsky 
- Zhabotinsky model, dn R-z=3.0 (0.02765- Dz Cm), Tp=20 (fp=0.05). Tp=20 
corresponds to 170sec. In this case the diode stops all pulses but the first one 
from the train in the forward direction, (curves 1 and 2). It works properly in the 


reverse direction (curves 3 and 4). 
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Figure 38. A snapshot showing the propagation of two consecutive pulses 
in a fragment of the cross junction structure. Brighter colors correspond to 
higher concentration of activator (u) The picture has been obtained for the 
FitzHugh - Nagumo model. The width of the passive stripes is 0.16 and the 
width of the active medium inside the channel equals to 1.02. The time shift 


between the pulses is £,=3.60 (f,=0.278). 
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the penetration depth for a pulse travelling along the stripes), but it is able 
to get through the perpendicular ”obstacle” (as the penetration depth in this 
direction is greater that the width of the passive stripes). However, the first 
pulse ”opens” the stripes which bound the channel for the following pulses. 
One can see that the next pulse ”leaks out” of the channel and excites the 
neighboring areas (cf. the second pulse in Fig. 38). Of course, for the 
FH-N model all the pulses in a train are able to cross the passive stripe 
perpendicular to the direction of motion. 

Similar computations have been carried out for the R-Z model with the 
width of the passive stripes of 3.291 (0.03034- Dzem) and the inner width 
of the active channel equal to 50.192 (0.4627. Dzem). They have shown 
that for signals of high frequency (eg. for 7,=20, corresponding to 170sec) the 
individual pulses propagate within the signal channel, but some of them may 
be stopped by the perpendicular passive stripe. Thus, the cross junctions 
described by both FH-N and R-Z types of dynamics are sensitive to the 


frequency of input signals. 


2.4.8 Discussion. 


This part of the work is concerned with the signals obtained after a regular 
train of pulses crosses a barrier of a passive medium. The most thorough 
calculations have been performed for the excitable medium described by the 
FH-N model, but qualitatively similar results may be obtained for the models 
of the BZ reaction (eg. the R-Z model and the Oregonator model considered 
here). 

Two facts are obvious: if the barrier is narrow - it is transparent to 
the pulses, if it is wide - it is impenetrable. However, it has been found that 


between these two limiting cases there is a range of barrier’s widths for which 
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it works as a transformer of signal frequency. It means that every second, 
third etc. pulse from the incoming signal is transmitted and all the others are 
stopped at the barrier. The number of transmitted pulses decreases with the 
barrier’s width. More complex examples of signal transformation have also 
been observed, like eg. a selection of two pulses out of every five arriving, 
shown in Fig. 20. Unfortunately, such interesting, complex behavior occurs 
in a narrow range of barrier’s width and it is not as robust as the simple 
division of the number of input pulses by two or three. 

The calculations have been performed for the FH-N model of an excitable 
system as well as for the R-Z and Oregonator models of the ferroin/cerium 
catalyzed BZ reaction. For all these models an interval of barrier’s width has 
been found in which the barrier works as a transformer of frequency of the 
input signal. The filtering ratio equal to 1/2 is dominant among the nontrivial 
filtering. It corresponds to elimination of every second pulse from the train 
by the barrier. In case of the FH-N model, this process may occur for barriers 
which are narrower or wider than the penetration depth for a single pulse 
(d.Fu-N=0.163). The ” scenario” of such elimination for barriers wider than 
d.FH-N in the FitzHugh - Nagumo model is shown in Fig. 17. In this case 
the first incident pulse ” dies” at the barrier, but the next one is transmitted, 
another one dies etc. Apparently a pulse stopped at the barrier increases the 
value of u inside the barrier for a short period of time, so the ” activator” is 
accumulated within the passive area. This helps the next incident pulse to 
get through the barrier. This mechanism is absent in the models of Belousov - 
Zhabotinsky reaction considered, where the frequency transforming has only 
been observed for barriers which are narrower than the penetration depth 
for a single pulse. In this case the stable elimination of every second pulse 


from a train means that the first incident pulse crosses the barrier, the next 
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one “dies” etc. For division by 3 - the first pulse gets through and then two 
pulses ”die” etc. (cf. Figs. 28 and 31).1% 

A family of transmitted signals has been found in which one signal out 
of n arriving gets through the barrier ((0)"~'(1) mode) . These signals, 
characterized by the firing number 1/n, contribute to the devil’s-staircase-like 
behavior of the firing number, expected on the basis of the ” naive” model of a 
perturbed excitable system (Section 1.2.2). However, less trivial patterns of 
the transmitted signal have also been observed. Another interesting family of 
transmitted signals, which can be described as (0)(1)"~!, corresponds to the 
case when only one signal is not transmitted out of every n arriving. Such 
behavior cannot be described by the ”naive” model, but comes out quite 
naturally, if one assumes that an excitable system which has not reached its 
stationary state may be re-excited by a perturbation which is strong enough 
(see discussion in [69]). Finally, many complex structures of the output signal 
have been observed. They follow the structure of the Farey tree [130, 131] 
(take for example points a - h in Fig. 22). 

From the point of view of numerical integration techniques, it is important 
to stress that the results concerning the complex transformations of chemical 
signals passing through a passive barrier are quite sensitive with respect to 


the time and space integration steps used. The exact position of plateaus in 


19]£ has been suggested [129], that the different behavior of the FH-N and R-Z models 
may be the result of different scales, within which the values of variables change in these 
models. Looking at Figs. 1 - 2 one may notice, that in case of the FH-N model the scale 
is the same for its both variables (u and v), while for the R-Z model the relative change 
of one of the variables (x) during the excitation cycle is several orders of magnitude larger 
than for the other variable (z). The author of this work has tried to adjust the parameters 
of the FH-N model to make its nullclines qualitatively similar to those of the R-Z model. 
Unfortunately, all attempts have been unsuccessful, so this hypothesis still needs to be 


verified. 
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the cross-sections like those presented in Figs. 22 - 24 depends on dt and 
dl, although the overall structure of the cross-section remains unchanged. In 
general, more exact results are obtained by means of more accurate algo- 
rithms (eg. adaptive grids), but even such algorithms seem to fail at the 
boundaries of plateaus (cf. Figs. 26 and 27) 

As a consequence of the frequency transforming phenomenon, it has been 
found that two dimensional reactors composed of active and passive regions, 
designed for direct processing of chemical signal, are sensitive to the input 
signal frequency. The operations performed on high frequency trains of pulses 
may differ from those performed on pulses which are well separated in time. 
This effect has been shown on two examples of such reactors: the chemical 
diode ([83, 85, 86, 102]) and the cross junction ([111, 112, 115]), but one 
should expect that the other signal processing devices (logical gates [86], 
memory device [87]) would also exhibit nontrivial dependence on the input 
signal frequency. Malfunction of the chemical diode, qualitatively identical 
with the results of calculations for the R-Z model presented in Fig. 36, has 
been observed experimentally by Y. Igarashi, K. Yoshikawa and J. Górecki 
at the Department of Physics, Kyoto University, Kyoto, Japan in spring 2002 
[132].20 

Transformation of chemical signal frequency on a passive barrier in the 
ferroin catalyzed BZ reaction has been recently reported by K. Suzuki, T. 
Yoshinobu and H. Iwasaki [70]. The diagram which relates the filtering ratio 
seen in their experiments with the barrier's width and the period of excita- 
tions (Figure 10 in [70]) is in a qualitative agreement with the results shown 


in Fig. 29. Moreover, individual examples of frequency transforming (in the 


20'Thanks to courtesy of the Authors, the film showing their experiment is included on 


the CD-ROM (diode_12.avi, folder Experiments) 
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ferroin catalyzed BZ system, too), have been observed experimentally by the 
author of this work. They correspond to several values of the filtering ratio 
(0, 1/2, 2/3, 1); please refer to Section 2.7 for details. 

The authors of [69] pointed out that the resonant patterns of transmitted 
waves may be important in biological systems. They postulated that narrow 
excitable gaps in an unexcitable tissue may be responsible for transforma- 
tion in the frequency of a biological signal. The results presented in this 
work show that such gaps are not necessary and the phenomenon may occur 
if some reagents responsible for signal propagation can diffuse through the 


unexcitable tissue. 
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2.5 Excitable pulses on a ring. 


Unidirectional waves on rings, described by FitzHugh-Nagumo equations, 
were used to study so called circumnutation of plants in their quest for a 
support [133]. A circular excitable medium with radial input and output 
channels can be used as a memory device [86, 87]. Thus, the features of ex- 
citable pulses propagating on a ring are important for applications of chemical 
systems to information processing, but also for the biologists. 

At first, a single travelling pulse on an interval [0,Z] with the periodic 
boundary conditions has been considered. It has been found [125] that the 
velocity, width and amplitude of the pulse in its stationary shape depend on 
L. The velocity and width grow remarkably with increasing L, reaching their 
maximum values, corresponding to propagation on an infinite interval. As an 
example, the results obtained for the FitzHugh - Nagumo model are presented 
in Figures 39 (velocity) and 40 (width). Here the width of the pulse is defined 
as the length of the interval in which the value of u is greater than 0.1 (cf. 
Figure 41). Changes in the amplitude of the pulse are not so remarkable. 
The amplitude increases together with Z by several percent only and quickly 
stabilizes. Real chemical systems are characterized by qualitatively the same 
dependencies between velocity or width of pulses and the space available for 
individual pulse. Consequently, several pulses initiated on a ring tend to 
arrange themselves in a symmetrical way, using the whole space available. 
This property of travelling waves was observed in real chemical experiments 
[89]. 

Once a pulse has passed through a part of the interval |0,Z] then a certain 
amount of time (called the recovery time) is needed before it would be able 
to propagate in this region again. If the pulse arrives earlier it cannot 


propagate and disappears. Thus, there exists a minimal length Lmin, below 
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Figure 39. The velocity of a single pulse in the FitzHugh-Nagumo model 
for various length L of a one dimensional system with periodic boundary 


conditions. Symbols mark the computational results. 
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Figure 40. The width of a single pulse in the FitzHugh-Nagumo model 
for various length L of a one dimensional system with periodic boundary 


conditions. Symbols mark the computational results. 
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Figure 41. The stationary shape of a pulse of u (” activator”) obtained for 
the FitzHugh - Nagumo model on an interval of length Z=20. The width of 
the pulse is defined as the length of the interval in which the value of u is 


greater than 0.1. 


CHAPTER 2. RESULTS. 119 


which the time, that the pulse needs to pass through the whole interval, is 
shorter than the recovery time of the medium. Consequently, for L< Lmin 
no stable travelling pulse can exist. After initiation the pulse on such a ring 
dies out. The critical length Lmin for one pulse has been found by numerical 
calculations. Then it has been confirmed that the critical length for two, 
three etc. pulses is just the corresponding multiplicity of Lmin. 

Using the implicit method with n=200, dt=5-10~? and t,,„=10000 it has 


been found that for the FitzHugh - Nagumo model 
Lmin, H-N = 3.529 (2.23) 


Similar calculations for the Rovinsky - Zhabotinsky model (n=200, dt=1-107%, 
Tmax =500) have shown that 


Dx 


Lmin.R- z X 5.79 (0.0534 - cm) (2.24) 


The value of Lmin can be easily transformed into the critical radius of a 
ring, supporting one stable pulse (called a ”1-pulse” ring for simplicity): 


| ee 
PE 225 
r % (2.25) 


The critical radii of rings, supporting two, three and more stable pulses are 
just corresponding multiplicities of Tmin. AS a consequence of Eqs. (2.23 - 
2.24) and Eq. (2.25) the critical radii of ”1-pulse” rings for the FitzHugh - 


Nagumo (Tmin.rH-N) and the Rovinsky - Zhabotinsky models (Tmin.R-z) are: 
Tmin,FH-N =~ 0.562 (2.26) 


Dx 


Tmin,R-z == 0.922 (0.00850 - cm) (2.27) 


Xo 
Two dimensional ring characterized by its outer radius R and inner radius 


r - as presented in Figure 42 - has been considered. The propagation of a 
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Dy 


Figure 42. The types of rings considered. A - a ring with R>rym and 
T>Tmin enables stable propagation of a pulse around it. B - a ring with 
R© Tmin and T<Tpmin Maker a rotating pulse change its shape during propaga- 
tion, because the pulse cannot constantly propagate along the inner boundary 


of the ring. 
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single pulse of excitation around the ring has been studied for various values 
of r and R. It is expected, that for r->rmin, R Tmin and R not much larger 
than r the pulse should just travels around the ring in a stable way, because 
in this case on the whole width of the ring the medium has enough time to 
relax before the excitable pulse reappears (cf. Fig. 42A). Consequently, one 
should use this kind of rings as memory devices, as proposed by I. Motoike 
and K. Yoshikawa [86, 87]. However, it is also expected, that for R>>ry;, and 
T<Tmin the pulse cannot constantly propagate along the inner boundary of 
the ring, because then the time of one cycle for a pulse travelling around the 
inner boundary is shorter than the refractory time of the medium (cf. Fig. 
42B). Thus, after one rotation the pulse has to ”leave” the neighborhood of 
the inner boundary, while it can still propagate closer to the outer boundary. 
Later, when the medium in the inner part gets relaxed again, the pulse may 
re-enter this area. As the result, a pulse rotating on such a ring becomes 
unstable. It is constantly changing its shape, moving towards and away from 
the inner boundary of the ring. 

The instabilities have been investigated for the FitzHugh - Nagumo model 
(Eqs. 2.1-2.2). The calculations have been carried out with the implicit 
method, adapted to use two dimensional polar coordinates. The grid con- 
sisting of 68 points in the radial direction and 512 points around the ring has 
been applied. In the calculations dt=5-10~° and tmar=130. No flux boundary 
conditions have been assumed on both boundaries of the ring. Pulses have 
been introduced to the ring using a prerecorded file, containing the station- 
ary shape of a pulse on a one dimensional interval with periodic boundary 
conditions. The mean values of u and v in the whole system (denoted by 
u and v) have been recorded as the functions of time and plotted in coor- 


dinates (u, v). Initial part of such trajectory, corresponding to transient 
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effects, has been ignored. The stable propagation of a pulse around the ring 
corresponds to a point on such a plot, as the values of u and y are constant 


1 Such stable propagation has been observed for R>T>Tmin.FH-N 


in time.” 
(in particular, R=3.0 or R=2.5 and r=0.57 have been used). On the other 
hand, when r<fmin,rg-n the plot becomes complicated, because the rotating 
pulse is constantly changing its shape, moving towards and away from the 
inner boundary of the ring. Here I have used R=3.0 or R=2.5 and r=0.53. 
An example of such trajectory (obtained for R=3.0 and r=0.53) is shown in 
Figure 43.2? For R=3.0 or R=2.5 and r=0.55 stable propagation of pulses 
around the ring has been observed. However, this is probably a transient 
effect, as the evolution of pulses on (two dimensional) rings has been studied 
up to tmac=130, while the value of Tmin.rH-N has been established in calcula- 
tions carried out up to tmarz=10000. Moreover, it has been observed that for 


R=2.0 pulses collapse after about two rotations, independently of the inner 


radius of the ring (r=0.53, 0.55 or 0.57). 


21To be very precise, this is in fact not a point, but rather a very small circle, because 
the rotating pulse in its stationary form is moving with respect to the grid used in the 


calculations; the diameter of this circle becomes smaller if a finer grid is used. 
22 An animation illustrating part of the evolution for which Fig. 43 has been obtained 


can be found on the CD-ROM (Fig_43a.htm, folder Figures). For comparison, stable 
propagation of a pulse around a ring with r=0.57 and R=3.0 is illustrated in Fig_43b.htm 
(folder Figures). 


CHAPTER 2. RESULTS. 123 


0.28 | | 
N 
NN 
INU 
UN 
| \ 


Average value of v 
SSS 


NZ R 


-0.20 -0.10 0.00 0.10 0.20 
Average value of u 


Figure 43. Average values of u and v (denoted by u and v in text) observed 
for a single pulse rotating on a ring with r=0.53 and R=3.0. The curve shown 
corresponds to te|30, 130]. The values of u and v are changing, because the 
pulse is constantly changing its shape while rotating. The result has been 


obtained for the FitzHugh - Nagumo model. 
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To sum up, the instabilities appear for TZTmin,FH-N, aS expected. On the 
other hand, they also seem to depend on R. The detailed mechanism leading 
to instabilities and the analytical condition for stable propagation of a single 


pulse on a ring are not known yet. This problem is still under investigation.”? 


23Even more interesting phenomena have been observed for a sphere with both ” poles” 
symmetrically cut off, so that the circumference along the ” equator” is much larger than 
Lmin, While the circumference along the ” parallels”, which remain after cutting off the 
poles, is smaller than Lmin. However, it has only been checked so far that there are 
several ” modes” of changing shape for a pulse rotating around such a sphere. The modes 
are characterized by different symmetry. The problem has not been investigated any 


further yet. 
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2.6 When two excitable pulses meet. 


Another very interesting type of complex evolution has been observed in the 
system composed of two active semiplanes separated by a passive stripe for 
the FitzHugh - Nagumo model [112]. It is presented in Figure 44, for which 
he same color convention as for Figs. 5 - 11 is used. The horizontal stripe 
of passive medium is transparent for a pulse propagating perpendicularly, 
but impenetrable for a pulse propagating parallel to it (cf. Fig. 44a). The 
first excitable pulse propagates along the passive stripe, to the right. The 
second pulse, initially located in the bottom part of the system, crosses the 
barrier and starts its propagation on the more relaxed side of the upper 
part, spreading to the right (following the first pulse) and upwards (cf. Fig. 
44b). On its way to the right it propagates along the passive stripe which 
has just been disturbed by the incident pulse. There is some activator left 
within the passive stripe, because the activator does not decompose there 
due to chemical reactions and it can only diffuse out of the stripe. The 
second pulse, propagating along the stripe to the right (cf. Fig. 44b), adds 
even more activator to the passive area (compare the ”whiskers” sticking 
out of the passive stripe, described in Section 2.3.1). Meanwhile, the active 
medium on the other side of the stripe becomes relaxed and ready for a new 
excitation. Because of that at some point the second pulse gets through 
the parallel barrier once again and starts spinning around (cf. Fig. 44c). 
While spinning, it is able to cross the passive barrier again and again (Figs. 
44d-f), due to the mechanism just described. Numerically the structure is 
stable. It persists in calculations lasting twenty times longer than a pulse 
needs to get through the whole square. No change in the position of the 
center of the pinwheel has been observed in the calculations carried out. 


This phenomenon has to be taken into account for cross junctions with very 
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Figure 44. The ”pinwheel” on a plane with a single stripe of diffusion field, 
obtained within the FitzHugh - Nagumo model. Gray areas correspond to 
the excitable field and the black stripe - to the diffusion one. Brighter colors 
correspond to higher values of u. The consecutive snapshots present the 
value of u at moments: a) t=4.0, b) t=7.0, c) t=10.0, d) t=12.0, e) t=14.0, 
f) t=200.0. 
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wide signal channels. 

The result illustrated in Fig. 44 has been obtained with the implicit 
method, on a 8x8 square covered with 400x400 grid. The width of the dif- 
fusion stripe here has been the same as for the cross junction described in 
Section 2.3.1 (dpry_n=0.16). The time integration step dt=0.005 has been 


used.?4 


24The animated version of Fig. 44 can be found on the CD-ROM (Fig_44a.htm, folder 


Figures). 
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2.7 The experiments on frequency transform- 
ing. 


The numerical results concerning the propagation of trains of pulses through 
a passive gap presented in Section 2.4 have been checked in experiments, 
run by the author of this work in the Institute of Physical Chemistry of the 
Wuerzburg University in Wuerzburg, Germany, in March 2002. The experi- 
ments were made in a system based on the ferroin catalyzed BZ reaction. 

Polysulfone membranes (PALL Corp., 47mm in diameter, 0.45um in 
width) were covered with ferroin (bathophenantroline, Fluka) exactly accord- 
ing to the recipe given by Lazar, Noszticzius, Forsterling and Nagy-Ungvarai 
in [89] (Sol3 and Sol4 in Section 2.1.2 of their article). At start, the BZ solu- 
tion providing the oscillatory regime for the reaction was prepared, accord- 
ing to [134] ([H2SO4]=1.25M, [NaBrO3|=0.07M, |CH2(COOH)2|=0.375M ). 
The rest of the experiment was carried out within the excitable regime, us- 
ing the BZ solution described by Suzuki, Yoshinobu and Iwasaki in [70] 
([H2504]=0.50M, |NaBrOs|=0.4M, |CH2(COOH)>2|=0.40M, [Na Br]=0.08 M). 
The results presented here were obtained for the excitable regime. 

In order to obtain two dimensional structures composed of active areas 
(with ferroin) and passive ones (without it) the elements of desired shape 
were first cut out of a ferroin covered membrane with scissors. The elements 
were then arranged on the adhesive tape and fixed. Alternatively, a piece of 
a ferroin covered membrane was first stick to an adhesive tape and later it 
was cut in two parts with a lancet or a razor blade (both techniques enabled 
to obtain very narrow passive stripes - of the order of 0.01mm - separating 
the active areas). Once built, the structures were fixed to the bottom of a 


Petri dish with the adhesive tape and immersed in the BZ solution. Usually 
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the experimental setup had the form of a ”comb”, providing several active 
(excitable) paths, along which pulses could propagate; on each path there 
was one passive barrier, ie. a slit between two pieces of a ferroin covered 
membrane.” Pulses of excitation were initiated by touching the membrane 
with a silver wire and they were destroyed or broken by pouring distilled 
water (water purified with Millipore Milli-Q filtering system was used for 
this purpose, as well as for preparing solutions; the resistivity of the water 
was 18.2MQ-cm). This way it was eventually possible to obtain a rotating 
spiral in a desired place of the system. The spiral was the source of a stable, 
regular train of pulses arriving at the passive barrier. The experiments were 
recorded with a CCD camera connected to a PC.”° 

The method of building the investigated systems did not allow to control 
precisely the width of passive gaps. However, once a device was constructed, 
it was possible to measure the distance between the active areas precisely 
under a microscope. Building a cross junction (as described in Section 2.3) 
was impossible with such a crude method of constructing the systems, so only 
the propagation of trains of pulses through a passive gap was investigated. 
The experimental observations are in qualitative agreement with the results 
presented in this work (Section 2.4): the frequency transforming effect has 
been observed.?” For example, filtering ratio f,/f,=2/3 was found for pulses 
of excitation arriving in about 60sec intervals at a passive barrier which 


was 0.016mm wide. It corresponds to elimination of every third of arriving 


25A picture showing one of such structures as an example can be found on the CD-ROM 


(Example.htm, folder Experiments) 
26]n such experiments the pulses are visible as blue stripes moving on a pink background. 


Blue and pink colors correspond to oxidized and reduced form of the catalyst (Fe(phen)3* 


and Fe(phen)3* respectively), compare [3] p. 183 or [4] p. 337. 
27Note that the concentrations of reagents used in these experiments were different from 


those assumed in the calculations for the R-Z model (Section 2.3.2). 
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pulses by the barrier.? The ” classical” division of the original frequency of 
incident pulses by 2 (f,/f,=1/2) was observed for pulses arriving every 46sec 
at a barrier 0.065mm in width.” Trivial examples of f,/f,=1 (all pulses 
get through the barrier) and f,/f,=0 (no pulse can cross the barrier) were 
observed for t,=26sec, d=0.016mm and t,=32sec, d=0.258mm respectively.” 

Moreover, the velocity of free propagation of a single pulse of excitation 
on a plain ferroin covered membrane was measured, giving 9.4mm/min. It is 
of the same order of magnitude as the value obtained from the calculations for 
the R-Z model (2.686mm/min for Dx=10~°cm?/sec, see Section 2.3.2), al- 
though the values of parameters used in computations (fixed concentrations, 


diffusion coefficients) were different from experimental ones. 


*8The recording of this experiment can be found on the CD-ROM (2-3.bat, folder Ez- 


periments) 
291 2.bat, folder Experiments on the CD-ROM. 
301 _1.bat and 0.bat respectively in folder Experiments on the CD-ROM. 


Chapter 3 


Possible applications of 
described devices to direct 


information processing. 


This section contains some ideas on potential future applications of reactors 
composed of active and passive areas to direct information processing. First, 
the advantages of simple structures of the cross junction and the passive 
barrier (Sections 2.3 and 2.4) are described. These two simple elements, 
together with the chemical diode, a logical OR gate and a memory cell (as 
proposed by I. Motoike and K. Yoshikawa - see Figures 2, 3 and 7 in [86] 
respectively) are combined together to build two more sophisticated devices. 
The first one is an "adding machine”, which is able to add two numbers 
given in binary representation. Second one is a counter of number of pulses 
that propagated in a signal channel. Both devices produce results in binary 
representation. Although some of the structures presented in this section 
may be hard to deploy, the author thinks that all of them can be constructed 


and applied in practice. 
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The cross junction with an extra channel, diagonal to the signal channels 
(cf. Figure 45A), may be easily used as a logical AND gate. One should 
concern paths 7, and i (Fig. 45A) as input channels and paths 01, 02 and 03 
as output channels. If and only if two pulses appear simultaneously in 7, and 
ło, then an output signal is produced in 03 (Fig. 5 is helpful to imagine what 
happens). If a single signal appears in 4 or k only - it propagates through 
the junction and goes out through o; or 02 respectively (and no signal is 
initiated in 03). Of course, there are no output signals at all if input signals 
do not appear in 4 nor ig. Thus, a logical AND gate has been constructed. 

The advantages of a simple passive barrier are even easier to show. First, 
a barrier of a fixed width d may be used as a filter of frequency of a chemical 
signal. If the frequency of the incoming signal f, is greater than a certain 
critical frequency feri (where ferit depends on d - cf. Fig. 29), then no pulse 
can cross the barrier. Thus, only signals with frequency f,<feri may be 
observed behind the barrier. The passive barrier may also be used to divide 
the number of pulses in the incoming train by a desired integer number. 
Knowing the frequency of pulses in the incoming train, it is enough to adjust 
the barrier’s width d appropriately (cf. Figs. 18, 19 and 29). The division by 
2 seems especially useful (it is used later in this section to build the counter of 
pulses), so it is very convenient that it appears in a wide range of parameters 
(d, f,). Last but not least, a ring ” equipped” with input and output channels 
may be used to regularize a train of pulses in which the time shift between 
consecutive pulses is random or unknown. First the whole train has to be 
” collected” on such a ring, then one should wait until the pulses arrange 
themselves symmetrically (this phenomenon is described in Section 2.5) and 
finally a regular train of pulses can be read from the ring. This procedure 


is limited by the ” capacity” of the ring used (so the longer trains are to be 
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Figure 45. Simple signal processing elements. A - the logical AND gate 
based on the cross junction. B - the logical OR gate. C - a memory cell (pulse 
storage ring). Thick black lines in the picture correspond to passive stripes 
of properly chosen width. Thin dotted lines in Figs. 45A and 45B resemble 
the symbols used to indicate the corresponding elements in the schemes of 


more advanced devices shown in Fig. 46. 
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regularized, the larger rings should be used). On the other hand, if a ”N- 
pulse” ring is used, then the frequency of the regularized train of pulses has 
to be within a certain (known) range. 

For the following I use two elements already built and tested. The logical 
OR gate (shown schematically in Fig. 45B) produces an output signal in o, if 
an input signal appears in i or e. The interference between the output signal 
and the input channels is eliminated with a diode-like structure present at 
the junction of 4, e and o,. This shape of the logical OR gate was proposed 
and confirmed to work in [86]. One more ”ready-to-wear” element which is 
needed is a memory cell (shown in Fig. 45C). It can be used to memorize 
and regularize a train of up to N pulses if a ”N-pulse” ring is applied. The 
pulses may later be read and/or erased. The structure presented in Fig. 45C 
was proposed and proven to work in model calculations as well as in real 
experiments in [87]. 

A scheme of a single segment of the adding machine is presented in Figure 
46A, in which the AND gates based on the cross junction (Fig. 45A) are 
indicated with squares, while circles stand for the OR gates (Fig. 45B). 
The structure shown in Fig. 46A is able to add the numbers” (”0” or ”1”) 
appearing in the input channels ap and bọ and memorize the result (in the 
form of a single pulse) in the memory cell Mp (thus, a” 1-pulse” ring is enough 
to serve as the cell). If a pulse appears in the input channel w (for simplicity 
let us write ag=1), it goes through the cross junction C4, then the OR gate 
R, and arrives at the cross junction (>. Here it may meet a pulse coming 
from bg. If it happens - both pulses annihilate within (> (so no output signal 
is produced in output channels 0; and 02 of (> - cf. Fig. 45A). Consequently 
no input signal appears on the OR gate R» and no signal (”0”) is memorized 


in My. However, an output signal is initiated in the output channel 03 of C2 
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(here the mechanism shown in Fig. 5 is used). This signal goes out of the 
device through the OR gate R3 and may by used as an extra input signal for 
the next segment of the adding machine (e1). The extra input channel eg is 
not really necessary in the first segment of the adding machine. However, its 
presence is crucial in all the other segments. An input signal appearing in 
this channel means that there was an ” overflow” when adding the numbers 
in the previous segment. Such ” overflow” signal should be added to a and 
b input signals in the next segment. The sum between a and e channels is 
performed in Cı in exactly the same way as the sum between dp and b in 
(>. The sum of a and e is then used as one of input signals for C2, where it 
is added to b. Using N segments like the one presented in Fig. 46A one is 


able to calculate the sum of two numbers given in binary representation: 


a= Get? a HE a +20! -aya (3.1) 


plus 


b = bo +2- by +4 bz +... +207. bya (3.2) 


gives 


M = Mo +2: Mı +4- M+... +207. My1 +2” - en (3.3) 


where all the a;, bi, M;€{0, 1}. 

A single segment of the counter of pulses propagating in a signal channel 
is presented in Fig. 46B. Before counting the original train of pulses should be 
collected on a ”N-pulse” storage ring (like the one shown in Fig. 45C). Thus, 
the counter is able to count up to N pulses. Once the train is regularized, it 
may be read from the ring and put as the input signal % to the first segment 
of the counter. The signal travels through the OR gate Ry and arrives at the 
passive barrier B. The width of the barrier should be properly adjusted to 


eliminate every second of the incoming pulses from the train (it is possible 
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Figure 46. Advanced direct information processing reactors. A - a single 
segment of the adding machine. B - a single segment of the counter. Squares 
symbolize logical AND gates based on the cross junction (cf. Fig. 45A), 
circles stand for the logical OR gates (cf. Fig. 45B). Lines mark the signal 
paths and arrows indicate the direction of propagation of pulses. The gray 


vertical interval in Fig. 46B marks the passive barrier. 
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because one is able to specify the range of frequencies of the regularized 
input signal). Let us also assume that the barrier’s width is set below the 
penetration depth for a single pulse, so that the first pulse from the incoming 
train is always able to cross the barrier. If there were K pulses in the initial 
train, then one would get [K/2]+1 pulses behind the barrier (where [p] stands 
for the ” floor” function, ie. the greatest integer number not greater than p). 
It is very important to check whether the initial number of pulses K was odd 
or even. This may be achieved by putting one more pulse of excitation in the 
channel e, just after the last pulse from the incoming train has arrived at the 
barrier B. The pulse from e is duplicated. First of its copies goes directly into 
the cross junction C4 and the second one goes through the OR gate Rı and 
arrives at the barrier B. If K was odd, then the second copy would ” die” at 
the barrier. The first copy would then arrive at Cı and go through it (in the 
vertical direction in Fig. 46B). Consequently a pulse would be memorized 
in the memory cell My, indicating that K was odd. Before memorizing, a 
copy of this signal is produced. This copy is directed to the cross junction 
C2, where it should be used to ”kill” the first of the pulses which previously 
passed through the barrier B (this is achievable if a proper difference in 
lengths of signal paths is provided). If K was even, then the second copy of 
e would cross the barrier and annihilate together with the first copy at the 
junction C4 (again, a proper difference of lengths of signal paths has to be 
provided). In this case no pulse would be memorized in the memory cell Mo, 
indicating that K was even and none of the pulses which previously passed 
through the barrier would be ” killed” at Cy. Anyway, the train of pulses 
going out of O» (in the horizontal direction), should be regularized (with a 
storage ring) and used as the input train for the next segment of the counter 


(i). The procedure should be iterated until all the pulses from the original 
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train (i) are eliminated. At his moment one gets the binary representation 


of K in the memory cells My, My, Mo etc. 
K = Mo+2- Mı +4: Mə +... (3.4) 


In case of the counter, we are not interested in the output signals which may 
appear in the diagonal channels of Ci or Ch (cf. Fig. 45A) when two pulses 
coming from perpendicular directions annihilate within the junction. Last 
but not least, it is worth mentioning that alternative structures of counters 
of pulses have already been proposed by J. Górecki, K. Yoshikawa and Y. 
Igarashi [135]. 

The adding machine and the counter would work properly if lengths of 
signal paths between their elements are chosen properly (so that the pulses 


would coincide in desired places at desired moments). 


Chapter 4 


Conclusions. 


The main purpose of this work was to study chemical reactors, composed 
of active and passive areas, from the point of view of their applications to 
direct processing of information in the form of chemical signals. 

Several types of two dimensional reactors have been studied. They include 
the cross junction (Section 2.3), the passive barrier (Section 2.4) and the 
ring (Section 2.5). All these structures have turned out to be useful in direct 
processing of chemical signals and thus the information that these signals 
may carry. The junction can be applied as a logical AND gate (coincidence 
detector of two pulses) or as a switch of the direction of propagation of a single 
pulse. It works for both FitzHugh - Nagumo and Rovinsky - Zhabotinsky 
models. 

Propagation of single excitable pulses as well as trains of such pulses 
through a passive barrier has been studied. It has been found that the 
barrier is able to systematically eliminate pulses from a train. Consequently, 
it can be used to divide the number of pulses in the train by a desired integer 
number. The mechanisms of elimination have been discovered to be different 


in the FitzHugh - Nagumo model and the Rovinsky - Zhabotinsky model, 
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in their forms and for the values of parameters used in this work (Eqs. (2.1 
- 2.4) and (2.5 - 2.8) respectively). However, it is not certain whether the 
observed differences come as the result of different parameters used, or is 
there a generic difference between the models. On the other hand, the fact 
that the elimination has been observed for three different models (FH-N, 
R-Z and the Oregonator) suggests that the phenomenon might be typical 
for all excitable systems, especially because it appears in the ”naive” model 
presented in Section 1.2.2. The phenomenon of systematical elimination of 
pulses from a train has turned out to be crucial for several signal processing 
devices (such as the chemical diode, logical gates etc.) previously proposed 
by other researchers, as well as for the cross junction proposed here (Section 
2.4.7). The results concerning stability of pulses on a ring may be important 
when such rings are used as memory cells. 

The numerical results concerning frequency transformations of regular 
trains of pulses by a passive barrier (Section 2.4) have been confirmed ex- 
perimentally. The results of experiments described in Section 2.7 are in 
excellent qualitative agreement with the numerical predictions. In particu- 
lar, frequency transformations corresponding to f,/f,=1, 2/3, 1/2 and 0 have 
been observed. These experimental observations are only some examples 
confirming the numerical results obtained. They imply that the Rovinsky 
- Zhabotinsky model used in this work describes pretty well systems based 
on the ferroin catalyzed BZ reaction. This conclusion is further supported 
by results of other experiments with chemical diodes (mentioned in Section 
2.4.7) and with passive barriers (eg. [70]). 

Two complex devices performing nontrivial operations, ie. adding two 
numbers or counting the pulses in a channel, have been proposed (Chapter 


3). Both of them utilize the basic signal processing elements, like the cross 
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junction and the passive barrier (discussed in this work) or the memory cell 
and the logical OR gate (known from the literature). The construction of the 
adding machine and the counter seems feasible, although it does not have to 
be easy. 

One final word of conclusion, coming from the last two paragraphs, is that 
once quite advanced information processing devices have been proposed and 
the model, within which these devices were ” born”, properly describes real 
chemical reactors (as the experiments suggest), we should expect to find such 
systems applied, for example in new generations of ” intelligent” materials or 


medicaments. Maybe in not so distant future. 
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Appendix A 


From the Field - Koros - Noyes 
mechanism of the BZ reaction 
to the Rovinsky - Zhabotinsky 


model. 


The contents of this part of the dissertation has been published as the sup- 
porting material for our publication in the Journal of Physical Chemistry A 
[115]. It is available at http://pubs.acs.org. 

When I prepared programs for numerical calculations I discovered that 
many papers which use the Rovinsky - Zhabotinsky model of the ferroin 
catalyzed Belousov - Zhabotinsky reaction contain simple mistakes in formu- 
las, which may significantly affect the results. Here I provide a step-by-step 
derivation of the dimensionless Eqs. (2.5-2.8) from the Field - Kórós - Noyes 
(FKN) reaction scheme and give the definitions of the scaling constants. 

Let us start from the modified Field - Körös - Noyes [42, 136] mechanism 


of the Belousov - Zhabotinsky reaction [44] completed by reaction (R13) (in 
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15); 

H* + HBrO3 + HBrO> = H BrO% + BrO- +H>O (R1) 

BrO,-+H* = HBrO; (R2) 

Fe(phen)3* + HBrO; = Fe(phen)3* + H BrO, (R3) 

2H BrO» > HOBr + HBrO3 (R4) 

Ht + Br + HBrO, — 2HOBr (R5) 

H+ + Br" + HOBr @ Br + HO (R6) 

H* + Br” + HBrO3 > HBrO; + HOBr (RT) 


Fe(phen)3* +CH Br(COOH); = Fe(phen)3* +H*+:CBr(COOH)> (R8) 


H.O +-CBr(COOH); > H+ + Br +-COH(COOH)2 (R9) 


HOBr + CH Br(COOH) © CBrz(COOH)> + H20 (R10) 
Br + CH Br(COOH),2 @ CBrą(COOH )> + H* + Br" (R11) 
CHBr(COOH) + H:O + CHOH(COOH))+H++Br- (R13) 


We use the notation as in [44]: A = [HBrOs3], B = |CHBr(COOH)>]|, 
C = [Fe(phen)3*] + [Fe(phen)3*], R = [-CBr(COOH),|, U = [H BrO}], 
X = [HBrO>], Y = [Br], Z = [Fe(phen)3*], ho is the Hammett acidity 
function [119, 120] and q is the stoichiometric factor. The kinetic equations 


introduced in [44, 45] read: 


X 
ox = —kyho AX +k_1U?+k3U(C—Z)—k_3X Z—2kahoX?—kghoXY +krhyAY 
(Al) 
oY 
OE = —kshoXY = krho AY + qkoR + kı3B (A2) 
OZ 
gy = BU(C — Z)- k-sXZ — kęBZ + k-shoR(C-Z) (A3) 
OU 


— = —2k1U? + 2kıho AX — ksU(C — Z) + k-3XZ (A4) 
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a = kgBZ — k_shoR(C — Z) — koR (45) 


where k4; are the rate constants of the corresponding reactions. 


These kinetic equations can be simplified. First, let us neglect the term 
kı: U in Eqs. (A1, A4) (as U is small [44]). Next, assuming that U and R 
are rapid variables comparing to X, Y and Z, from Eq. (A4) we get: 


_ 2hqhoAX + k_gXZ 


ka(C- Z) sy 


and from Eq. (A5): 


kg BZ 
R= AT 
k_gho(C — Z) + kg s 


In the denominator of the last formula we can neglect kg, because koKk_ghg(C- 


Z). All these approximations are based on experimental observations, quoted 
in [44]. 
Substituting (A6, A7) to (Al - A3) one obtains: 


X 
ox = ki ho AX — 2kgho.X? — kshoXY + krhoAY (A8) 
OY KgBZ 
= — kshoXY — krhoAY B A 
a ne-a eS A AE Sz 
ƏZ K;BZ 
— = ih AX —— => A10 
g ah ho(C — Z) van 
where 
ea (A11) 
k-s 


Up to this point the procedure presented here is identical as in [44]. Basing 
on Eqs. (A8 - A10) we derive the dimensionless reaction - diffusion equations 


for active and passive fields. 
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For the active field we eliminate Y as a fast variable using Eq. (A9) (see 


[44]) and get: 


K;BZ 1 
Y = (g———— + ks B)- A12 
(ocz PCR RXTEKA PA 
which can be substituted to Eqs. (A8) and (A10). It follows that: 
Ox k5X —k7A KgBZ 
— = ky hg AX — Zk4hy X” — : + kiz B Al 
Ot 170 aho EETA (CoD) 13B) (A13) 
OZ K;sBZ 
— = ŻkyhoAX — ——— Al4 
at ati BCE one 


In the passive field we have Z=C=0 and so R=0. This reduces Eqs. (A8 - 


A9) to: 
0X 


gy = —ŻkahoX? — kshoXY + krhoAY (A15) 
oY 
Sp = KsłoXY — krhoAY + kisB (A16) 


Once again the fast variable Y is eliminated. From (A16) we receive: 


s ki3B 
— kshoX + krhoA 


(A17) 


and after substitution to (A15) we end up with the following kinetics for the 


passive field: 


OX k5 X = kz A 

ESEE RY RE p Al 

a kaho EN KA (A18) 
Z=0 (A19) 


In our model the activator X diffuses, while the catalyst Z is immobilized. 
Thus, we add the diffusion of X to the kinetic equations (A13-A14) for the 


active and (A18 - A19) for the passive regions and get: 


X 
ox == ky ho AX = Zkąhy X? ze ( 


kX = kA) KsBZ 
kgX +kyA’  "ho(C-Z) 


+ k13B) + DxV2X 
(A20) 
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OZ KgBZ 
— = 2kj hp AX — ——~ A21 
Ot R ho(C — Z) aż 
for the active field and 
OX k5X —k7A 
= = =2kuhyX” — ky3B + DxV2X A22 
= kaho XTRA” xV; (422) 
Z=0 (A23) 
for the passive one. 
If one introduces dimensionless variables z, z, 7 and p, so that 
kA 
Xess A24 
ok, (A24) 
Z=Cz (A25) 
k4C 
= A26 
kŁAŻhy | ee) 
Ome | 
= gp eS ATs a A27 
r by A x':p ( ) 
and coefficients: 
ki KB 
= A28 
= RÆ 2 
2k4k13B 
z RZ A29 
2kąkz 
= A30 
R (430) 
kA 
eee A31 
kC (431) 
then Eqs. (A20 - A21) and (A22 - A23) turn into: 
OF Soa Oly he (432) 
ðr e Bee a e 
Oz z 
BE eto M A33 
OT A 31 -z ( ) 
for the active medium and 
1 = 
ea a, (434) 
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z=0 (A35) 


for the passive medium. Eqs. (A32 - A33) and (A34 - A35) are identical 
with Eqs. (2.5-2.6) and (2.7-2.8) in Section 2.1.2, respectively. 

Now let me point out the differences between my results and the other 
papers. The scaling of time (Eq. (A27)) and definition of 8 (Eq. (A29)) differ 
from those given in [44, 45] by the factor ho in the denominator, although 
the remaining scaling constants (Eqs. (A24), (A25), (A28), (A30), (A31)) 
and the reaction - diffusion equations (A32 - A33) are identical with those 
given in [44, 45]. 

The Rovinsky - Zhabotinsky model was also used by Kusumi et. al. in 
[84]. Their kinetic equations for the active field are identical with (A32 - 
A33), however the definition of u they give differs from (A30) by the factor 
of 2 and the scaling of p is different from (A27) by the factor ha! * in the 
denominator. The passive medium in [84] is described by kinetic equations 
simpler than (A34 - A35). The simpler equations were however not sufficient 
for my study, because for the values of parameters considered in this work 
they lead to excitations of active areas at the boundary of active and passive 


media. 


Appendix B 


Integration schemes used. 


Let us consider a reaction - diffusion equation in one dimension in the most 
general form (B1): 
— = D— +r(c) (B1) 
z 


where c is the concentration of a chemical reagent, t and x are the time 
and space variables respectively, D stands for the diffusion coefficient of the 
reagent and r(c) is the kinetic term, describing the chemical reaction. 

Let us assume that the values of c are discretized on a grid. Let us denote 
the current value of c at the i-th grid point (4=0,1,...,n,n+1) by cj” (known) 
and the corresponding value after time dt by ct (unknown). Let us assume 
that the grid is uniform and the distance between consecutive grid points is 
dz. 


In general cj” can be calculated as: 


cj =d + Ac; + dt Bo] (B2) 
i “i 1 Ox2"" 


where Ac; corresponds to the change of c at the i-th grid point due to the 
chemical reaction. 


In the Euler explicit method the calculations of c;” are fully based on the 
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current state of the system c;”. Consequently: 


8?c 


f— L +dt-r(d)+dt- D i B3 
ci = e + dt- r(e?) + dt- DES] (B3) 

where: 
OFE, 644 Poj 234 (B4) 


ox? i= (dz)? 
In the Crank - Nicolson method the diffusion term is approximated as: 


2 0 0 0 t t t 
O*c oi Gag PGA 258 VGT ECE iS, 


o. (dz)? | (dz)? 


| (B5) 


This makes c; occur on both sides of Eq. (B2) and shows that the methods 
based on the Crank - Nicolson discretization of the Laplace operator are 
implicit methods. 


Equation (B2) can be solved with respect to ct giving: 


CG; = Ai- Gi, +B; 1=0,1,...,n (B6) 
where Dg 
Ą = CZ (B7) 
1- (dz): (Ai-1 - 2) 


ae? (Bia + cja + Ch — 27) + Ac + ef 


1- PEA -2) 


B; = 


(B8) 


Additional constraints are imposed on the coefficients A;, B; by the boundary 
conditions used. 

In case of no flux boundary conditions we have Aọ=1 and Bo=0 (from 
Eq. (B6), because co’=c,'). Consequently, all the A;, B; up to An, Bn can 
be calculated using Eqs. (B7 - B8). At the other boundary of the system, 
according to Eq. (B6) we have: 


c = Ån Chai + Ba (B9) 
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but due to no flux boundary conditions c,'=c,+1°, SO: 


t Bn 
Cn+1 = 1— A,n, 


(B10) 


Having cn41*, the values of c;* for i=n,n-1,...,1,0 can be calculated by back- 
substitution, using Eq. (B6). 
In case of periodic boundary conditions, one has to solve a more compli- 


cated set of equations: 


ch = Ao: cq + Bo (Bila) 

dh = A h+ By (B11b) 
ch = Ana: É + Bri (Blic) 
G = ch (Bild) 

ApS Arme = Ana (B1le) 


Equations (Blld - Blle) come directly as the consequence of the periodic 
boundary conditions (B11d) and the fact that in this case no grid point on 
the interval is distinguished (Blle). 

In Section 2.4.1 a modification of the implicit method just described is 
mentioned. Instead of a uniform (equally spaced) grid, an adaptive grid 
is used. The adaptive grid is characterized by a smaller distance between 
consecutive grid points within and near the passive barrier and longer one 
far from the barrier. Let us denote the space step to the left from the i-th 
grid point (ie. the distance between nodes no. i-1 and i) by dz;_; and the 
space step to the right from the i-th grid point (ie. the distance between 
nodes i and i+1) by dz;. 


APPENDIX B. INTEGRATION SCHEMES USED. 162 


Then the Laplace operator is given by Equation (B12): 


2 0 0 t t 
O c 1 2 Ĉi—1 | Ĉi+1 d£zi-ı + dx; OW C; 1 Cig, dti 1 + dz; t 


Jr = 2 dr; ı + dx; dx;_4 | dz; dx;_\dx; RZA ie | dz; dz; dz; c] 
(B12) 


Equation (B6) holds and consequently the coefficients A;, B; have more com- 


plicated forms: 


Dd 1 
— dzi—1+dzi dx; 
Ai = | — — Ddt (a= STU bean ) (B13) 
dx;—,+dx; d£zi—1 dx; d£zi—ı 
Ddt ($u | Sae | "a LAG + 2 
; : - I 7 I | I Ci T Ci 
B = dzi—1ı+dzi dx; dtp 4 dzi—ı (B14) 
nrz | — Pát c= magii ee ) 
dx;—,+dx; dzi dx; dzi—ı 


The kinetic term Ac;, present in Eqs. (B8, B14), can be calculated in 
many ways. For the purpose of this work, the well established Runge - 
Kutta method of the fourth order has been used. For computations in two 
dimensions (eg. the cross junction - Section 2.3, the rings - Section 2.5 and the 
pinwheels - Section 2.6) the implicit method presented here has been applied 
to individual rows of the two dimensional grid and then the procedure has 
been iterated to obtain self-consistency between the rows. 

Ddt 


The explicit Euler method becomes unstable when 224 > 1 


CHE 5, So for a 


fine grid a very short time integration step is required, whereas the implicit 


method is stable for all values of dz and dt. 
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Contents of the CD-ROM. 


The CD-ROM contains supporting material for the PhD thesis. It has the 


following structure: 
readme.htm - table of contents of the CD-ROM. 
thesis.pdf - the full version of the PhD thesis in PDF format. 


\Experiments - films showing selected frequency transforming modes 
observed in experiments with the ferroin catalyzed Belousov - Zhabotin- 
sky reaction. Double-click the .BAT files to see the corresponding 
films." The file diode_12.avi shows malfunction of the chemical diode, 
observed experimentally by Yasuhiro Igarashi, Kenichi Yoshikawa and 
Jerzy Górecki at the Department of Physics, Kyoto University, Kyoto, 
Japan in spring 2002.? 


\Figures - animated extensions of selected figures from the thesis and 


some additional animations. 


1Some computers with Windows NT 4.0 Workstation fail to display the films. 
2Courtesy of the Authors of the experiment. 
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\Extra 


\Acrobat - Acrobat Reader 4.05 installation program for Mi- 
crosoft Windows 95/98/NT/2000/XP. Acrobat Reader is a tool 
for viewing PDF files. To install it, double-click the rs405eng.exe 
file. 


\Cross_junction - a complete web presentation about the cross 
junction within the FitzHugh - Nagumo model. To start browsing, 
double-click the !START.htm file. 


* All documents on the CD-ROM are optimized for viewing with Mi- 
crosoft Internet Explorer 4.0 or later, in 800x600 (256 colors) full screen 
display mode. In order to switch to the full screen mode start Internet Ex- 
plorer and select View >>> Full Screen or just press F11. To switch off the 


full screen mode press F11 again. 
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